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In  the  point-mass  algorithm  which  would  compromise  the  accuracy  in  the  deter¬ 
mination  of  the  oceanic  geoid,  attainable  with  the  SEASAT  altimetry  or  with  a 
similar  high-precision  observational  system. 

A  mathematical  framework  Is  provided  in  order  to  model  observations  such 
as  the  horizontal  (north  and  east)  and  vertical  gradients  of  gravity  or  gravity 
disturbance.  This  development  is  based  on  the  tensor  approach  to  theoretical 
geodesy  and  is  mostly  concerned  with  the  second-order  derivatives .along  local 
Cartesian  axes,  of  a  general  scalar  function  of  position;  the  derivatives  are 
expressed  in  spheroidal  and,  for  a  variety  of  practical  applications,  in  spheri¬ 
cal  coordinates.  The  scalar  function  of  position  is  specialized  to  represent 
the  (total)  earth's  potential,  the  standard  (or  normal)  potential  and,  especial¬ 
ly,  the  disturbing  potential.  ^ 

This  specialization  allows,  among  other  things, to  model  the  gradients  of 
gravity  disturbance  directly  by  a  suitable  set  of  parameters,  such  as  spheri¬ 
cal  harmonic  potential  coefficients  and  point  mass  magnitudes.  The  pertinent 
adjustment  algorithms  developed  in  terms  of  both  these  kinds  of  model  para¬ 
meters  are  presented  in  detail.  In  addition,  similar  algorithms  are  recapitu¬ 
lated  for  geoid  undulations  and  gravity  anomalies,  and  developed  for  the  de¬ 
flections  (north  and  east)  of  the  vertical.  The  adjustment  algorithms  in  terms 
of  the  potential  coefficients  contain  minimal  approximations,  while  those  in 
terms  of  the  point  masses  contain  the  familiar  spherical  approximation. 

Two  further  topics  discussed  in  this  report  are  related  to  the  sea  sur¬ 
face  topography.  The  first  deals  with  the  average  influence  of  the  astronomi¬ 
cal  tidal  forces  on  the  flattening  of  an  idealized  mean  sea  level,  and  the 
second  is  concerned  with  the  development  of  a  model  in  which  the  presence  of  sea 
surface  slopes  with  respect  to  an  equipotential  surface  is  taken  into  account. 

In  a  separate  development,  a  comparative  analysis  is  carried  out  for  the  node¬ 
point  approach  and  the  point-mass  approach,  focusing  on  the  deterministic  par¬ 
tial  derivatives  entering  the  observation  equations. 
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1.  OUTLINE  OF  THE  RESEARCH  TOPICS 


In  the  recent  past,  research  has  been  conducted  aimed  at  im¬ 
proving  the  knowledge  of  the  earth's  gravity  field  and  its  fundamental 
surface,  the  geoid.  The  most  extensively  treated  observational  mode  has 
been  satellite  altimetry.  The  adjustment  model  utilized  in  this  task  has 
been  the  short-arc  model,  described  and  refined  in  several  reports  and 
papers,  e.g.  in  [Brown,  1973]  and  [Blaha,  1975,1977,1979].  The  body  of 
altimeter  data  has  been  initially  obtained  via  the  GEOS-3  satellite.  How¬ 
ever,  several  aspects  of  the  SEASAT  observational  system  have  been  investi¬ 
gated  along  with  the  GEOS-3  system  in  view  of  improved  determinations  of 
the  oceanic  geoid.  For  example,  the  SEASAT  altimeter  has  been  considered 
in  conjunction  with  the  so-called  precise  ephemeris  (having  the  standard 
error  in  position  of  approximately  Zm)  and  the  recommendation  has  been 
reached  that  the  time  interval  for  each  satellite  arc  should  not  signifi¬ 
cantly  exceed  7  minutes,  in  order  to  eliminate  aliasing  effects  that  could 
become  noticeable  at  points  far  from  the  mid-arc  epoch. 

The  geoidal  surface  or  its  parts  have  been  adjusted  in  two  steps. 
The  first  adjustment  has  proceeded  in  terms  of  a  truncated  set  of  spherical 
harmonic  potential  coefficients  involving  global  data  sets  and  resulting  in 
a  fairly  smooth  global  geoid.  The  second  adjustment,  based  on  regional  al¬ 
timetry  residuals  from  the  first  adjustment,  has  had  point  mass  magnitudes 
as  parameters  and  its  main  purpose  has  been  a  more  detailed  (but  localized) 
determination  of  the  oceanic  geoid.  The  location  of  the  point  masses,  in¬ 
cluding  their  depth  underneath  the  earth's  surface  considered  as  a  sphere. 


has  been  stipulated  beforehand.  A  question  has  arisen  whether  the  point 
mass  parameters  could  serve  for  both  the  global  and  local  geoid  determina¬ 
tions.  However,  the  discussion  in  Chapter  2  shows  that  due  to  accuracy  con¬ 
siderations  such  an  approach  is  not  desirable.  In  particular,  the  geoid  de¬ 
termination  could  be  affected  by  errors  of  several  decimeters  due  to  the 
spherical  approximation  in  the  point  mass  algorithm  alone.  On  the  other 
hand,  if  the  point  mass  adjustment  is  based  on  a  higher-order  surface  as 
obtained  through  a  global  spherical  harmonic  adjustment  this  approximation 
will  not  compromise  the  accuracy  of  the  results. 

Originally,  the  satellite  altimetry  has  been  complemented  by  mean 
gravity  anomalies  which  have  prevented  the  system  of  normal  equations  from 
being  ill-conditioned  due  to  the  gaps  in  altimeter  data  (especially  over 
the  continents).  Eventually,  however,  need  has  been  felt  for  the  modeling 
of  additional  kinds  of  observations.  In  order  to  satisfy  such  requirements. 
Chapter  3  provides  the  mathematical  framework  for  several  observational  modes, 
notably  for  horizontal  and  vertical  gradients  of  gravity  or  gravity  distur¬ 
bance.  It  is  mostly  concerned  with  the  development  and  application  of 
second-order  derivatives,  along  local  Cartesian  axes,  of  a  scalar  function  of 
position,  expressed  in  spherical  coordinates.  This  development  is  based  on 
the  tensor  approach  to  theoretical  geodesy,  introduced  by  Marussi  [l95l] 
and  further  elaborated  by  Hotine  [1969],  which  can  lead  to  significantly 
shorter  demonstrations  when  compared  to  conventional  approaches. 

The  second-order  derivatives  of  Chapter  3  are  formulated  in  quite 
general  terms  and  are  applicable  to  any  scalar  function  of  position  (denoted 
F).  If  F  becomes  W  (the  earth's  potential)  the  six  distinct  second-order 
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derivatives  --  which  include  one  vertical  and  two  horizontal  gradients  of 
gravity  —  relate  the  symmetric  Marussi  tensor  to  the  curvature  parameters 
of  the  field.  The  second-order  derivatives  are  developed  also  in  spheroidal 
coordinates.  This  serves,  besides  its  own  theoretical  interest,  for  a  veri- 
ficiation  of  the  desired  formulas  in  spherical  coordinates  by  confirming 
them  indirectly,  through  a  coordinate  transformation  involving  differential 
quantities  and  relationships. 

The  general  formulas  for  the  second-order  derivatives  of  F  are 
specialized  to  yield  the  second-order  derivatives  of  U  (standard  potential, 
often  called  normal  potential)  and  of  T  (disturbing  potential),  which  allows 
the  latter  to  be  modeled  by  a  suitable  set  of  parameters,  for  example,  by 
spherical  harmonic  potential  coefficients,  chosen  localized  parameters  such 
as  point  mass  magnitudes,  etc.  The  second-order  derivatives  of  T  where  the 
property  AT = 0  is  explicitly  incorporated  are  also  given.  According  to  the 
required  precision,  the  spherical  approximation  may  or  may  not  be  desirable; 
both  kinds  of  results  are  presented.  The  derived  formulas  can  be  used  for 
the  modeling  of  the  second-order  derivatives  of  W  or  T  at  ground  level  as 
well  as  at  moderate  or  high  altitudes.  They  can  be  further  applied  in  a  ro¬ 
tating  or  a  nonrotating  field. 

In  Chapter  4,  adjustment  algorithms  are  developed  for  five  ob¬ 
servational  modes  including:  geoid  undulations,  gravity  anomalies,  deflec¬ 
tions  (north  and  east)  of  the  vertical,  horizontal  gradients  (north  and  east) 
of  gravity  disturbance  and  vertical  gradients  of  gravity  disturbance.  Each 
model  is  presented  in  the  spherical  harmonic  (S.H.)  formulation  and  in  the 
point  mass  (P.M.)  formulation.  The  former  does  not  contain  the  familiar 
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spherical  approximation  of  the  earth,  while  the  latter  does.  The  procedure 
for  computing  predicted  values  of  the  above  observation  modes,  as  well  as 
their  variances,  is  a  part  of  this  development.  Also  included  is  the  practi¬ 
cal  discussion  regarding  the  cut-off  distance  in  the  P.M.  adjustment,  i.e., 
the  size  of  a  spherical  cap  centered  at  an  observation  point  beyond  which 
the  P.M.  parameters  are  ignored  during  the  formation  of  an  observation  equa¬ 
tion. 

Of  the  above  five  models,  the  first  two  (i.e.,  for  geoid  undula¬ 
tions  and  gravity  anomalies)  have  been  developed  earlier.  However,  they 
are  recapitulated  in  a  concise  manner,  with  the  cut-off  considerations ’add¬ 
ed  on,  and  their  format  conforms  to  that  of  the  other  three  models.  Most 
of  the  new  development  is  based  on  the  formulas  of  Chapter  3.  In  this  re¬ 
spect  Chapter  4  is  the  only  chapter  built  upon  another.  Most  of  the  chapters 
in  this  report  as  well  as  the  appendix  are  presented  essentially  in  a  self- 
contained  manner  and  can  be  read  independently.  A  good  example  of  this  for¬ 
mat  is  Chapter  3  whose  shortened  version  is  being  presented  separately  as 
[B1 aha,  1980], 

Chapter  5  contains  two  topics  related  to  the  sea  surface  topography. 
The  first  deals  with  the  average  influence  of  the  astronomical  tidal  forces 
on  the  flattening  of  an  idealized  mean  sea  level,  and  the  second  is  concern¬ 
ed  with  the  development  of  a  model  in  which  the  presence  of  sea  surface  slopes 
with  respect  to  an  equipotential  surface  is  taken  into  account.  Such  slopes 
are  known  to  exist  in  different  parts  of  the  world  where  they  can  reach  a 
few  decimeters.  Considering  the  precision  of  the  SEASAT  observational  system 
or  other  such  high-precision  systems,  these  features  are  detectable.  They  can 
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be  modeled  and  eventually  incorporated  into  the  adjustment  process.  By  re¬ 
peating  this  process  at  different  time  epochs,  the  time  history  of  sea  sur¬ 
face  slopes  in  an  area  such  as  the  Equatorial  Pacific  Ocean  can  be  studied. 


The  appendix  of  this  report  is  mostly  concerned  with  comparing 
the  node-point  approach  and  the  point-mass  approach.  The  bulk  of  such  a 
comparison  —  which  could  include  a  multi-quadric  model  or  other  models 
used  to  describe  the  local  geoid  —  is  the  study  of  the  "covariance"  func¬ 
tion  giving  rise  to  the  deterministic  partial  derivatives  forming  the  ma¬ 
trix  of  observation  equations.  The  analysis  indicates  that  the  main  dif¬ 
ference  between  the  above  two  approaches  consists  in  the  somewhat  different 
shapes  of  their  respective  "covariance"  functions.  If  the  depth  of  the 
point  masses  corresponds  to  the  complexity  of  the  reference  surface  form¬ 
ing  the  basis  of  the  node-point  approach,  the  similarity  between  the  de¬ 
terministic  partial  derivatives  in  the  two  approaches  can  be  remarkably 
good. 
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2.  SPHERICAL  APPROXIMATION  AS  A  LIMITING  FACTOR 
TO  USING  POINT  MASSES 
IN  A  GLOBAL  GEOIDAL  REPRESENTATION 


A  point  mass  model,  in  which  the  adjustable  parameters  have  been 
the  point  mass  magnitudes  while  the  location  of  the  point  masses  has  been 
stipulated  beforehand,  lends  itself  to  an  efficient  treatment  only  if  the 
spherical  approximation  is  part  of  the  model.  This  has  been,  in  fact,  the 
approach  followed  in  previous  applications  which  have  led  to  a  detailed 
geoidal  representation  over  a  limited  area.  In  such  applications,  alti¬ 
meter  data  collected  by  GEOS-3  satellite  have  been  utilized,  together  with 
the  so-called  broadcast  ephemeris  (having  a  positional  standard  error  of 
about  20  meters).  A  question  emerges  when  a  global  geoidal  resolution  -- 
as  opposed  to  a  localized  resolution  --  is  contemplated  in  terms  of  the 
point  mass  parameters  alone:  Would  the  potential  accuracy  of  such  a  model 
be  acceptable  in  view  of  the  SEASAT  (or  a  similar)  observational  system 
which  is  far  superior  to  the  GEOS-3  system  (in  both  altimeter  and  ephemeris 
accuracy)? 

We  shall  first  consider  a  few  aspects  in  the  formulation  of  a 
global  point  mass  model  and  then  discuss  the  applicability  of  the  model  in 
terms  of  future  accuracy  requirements  approaching  one  decimeter  in  geoidal 
determination.  At  the  outset,  such  a  model  should  be  subject  to  the  global 
constraints  defining  the  coordinate  system.  In  spherical  harmonic  applica¬ 
tions,  the  coordinate  system  has  been  usually  defined  by  stipulating  the 
values  of  the  following  spherical  harmonic  potential  coefficients: 
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(2.1) 


cio=0»  Cir°*  Sji-O.  C2i=0,  S2i=0,  S22  =  constant. 

The  first  five  relations  correspond  to  the  "forbidden  harmonics"  and  the 
value  of  the  last  one  depends  on  the  orientation  of  the  x,y  axes  with  re¬ 
spect  to  the  principle  axes  of  inertia. 

The  following  expressions  are  given,  in  a  similar  form,  on  page 
61  of  [Heiskanen  and  Moritz,  1967]: 


iQ0  =  /  d(kM)  *  kM  ; 

(2.2a) 

i10  =  /zd(kM),  Ajj  *  /  xd(kM), 

B1X  -  /  yd(kM)  ,  1 

(2.2b) 

i21  =  /xzd(kM),  B21  - /yz  d(kM), 

B22  =  H  /xyd(kM)  ;  J 

i20  =  h  f  (2z2  -  x2  -  y2)  d(kM)  ; 

(2.2c) 

i22  *  *  /  0<2-y2)  d(kM)  • 

(2. 2d) 

Here  the  integral  extends  over  the  whole  earth,  k  is  the  gravitational  con 
stant  and  M  is  the  earth's  mass.  Further  we  have  (see  Section  2-5  of  the 
same  reference): 


A™ 

Anm 


’  <kM>  a"  c™ 
’  (kM>  a"  snm 


» 


9 


where  a  is  the  earth's  equatorial  radius. 
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Considering  a  discrete  distribution  of  mass,  the  formulas 
(2.2b)  -  (2. 2d)  yield 


10  =  (kM)  a 

C10 

1^0  •r-i 

II 

zj(kH)j  • 

u  =  (kM)  a 

C11 

n 

*imj  ■ 

U  5  (kM)  a 

S11 

-  I 

j 

yj(kM)j  . 

(2.2b') 


A21  s  (kM)  a2  C21 
B21  =  (kM)  a2  S21 
B22  =  (kM)  a2  S22 

A20  H  ^  a*  C20 
An  =  (kM)  a2  CZ2 


-  I 

j 

-  I 

j 

=  hi 

D 

*  hi 

0 

*  hi 

j 


yjZ.fkMlj  , 


XjyjtkHlj  ; 


J 


(22^-kj-y])  (kM)3  i 
(xj-yj)  ( kM) j  . 


(2.2c') 
(2. 2d') 


From  (2.2b*), 

I  z<(kM). =  0  , 
j  J  J 

?  xjzj(kM)r° 

3 


the  relations  (2.1)  now  read 


‘(2.3) 

’  I  yj2j(kM)j  =  0  ’  t  xjyj(kM)j  *  2(kM)  a*  S22* 

<5  J 
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Further, 

l  (2zj'Xj'yp  *kM*j  *  2{kM^  a*  C20  ;  (2.4) 

0 

I  (x’-yj)  (kM)j  =  4(kM)  a2  C22  .  (2.5) 

J 

However,  it  is  more  advantageous  to  deal  with  the  disturbing 
potential  (T)  than  with  the  actual,  "full"  potential  and  the  related  quan¬ 
tities,  since  by  far  the  largest  point  mass  magnitudes  would  be  "wasted" 
in  reproducing  the  known  normal  field.  Since  all  the  potential  coeffi¬ 
cients  except  C2Q  ,  C4g  ,  etc.,  are  zero  in  the  normal  field,  the  equations 
(2.3)  defining  the  coordinate  system  apply  exactly  as  they  stand  also  in 
this  new  situation.  In  the  adjustment  terminology,  they  represent  con¬ 
straints.  If,  in  a  point  mass  adjustment,  an  equivalent  of  the  ellipsoidal 
coefficient  C2Q  should  be  enforced,  the  relation  (2.4)  would  represent 
another  constraint  whose  right-hand  side,  however,  should  be  set  to  zero 
(it  would  now  represent  the  departure  from  the  value  corresponding  to  the 
ellipsoidal  C2q).  If  an  equivalent  of  C22  should  be  held  at  a  predetermined 
value,  the  relation  (2.5)  would  become  a  constraint  exactly  as  it  stands. 

Besides  the  six  constraints  presented  in  (2.3)  —  and  perhaps  two 
or  more  optional  constraints  similar  to  (2.4),  (2.5),  etc.  —  two  additional 
constraints  could  be  utilized  In  order  to  match  the  constraints  often  en¬ 
countered  in  adjustments  of  spherical  harmonic  potential  coefficients.  In 
particular,  the  mass  of  the  reference  ellipsoid,  taken  as  the  mean  earth 
ellipsoid,  and  its  potential  are  usually  preserved  throughout  the  adjustment. 
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This  corresponds  to  preserving  the  following  relationship: 


Jd(kM)  =  0  , 

/  T  do  =  0  , 

where  the  first  integral  extends  over  all  of  the  "anomalous  masses"  and  the 
second  integral  extends  over  the  whole  unit  sphere.  In  terms  of  the  point 
mass  parameters,  these  two  equations  represent  the  constraints  of  the  fol¬ 
lowing  form: 

I  (kM)j  =  0  , 
j 

l  ?  *1  * 0  ■ 

1  J 

where  l. -  is  the  distance  between  the  j-th  point  mass  and  the  i-th  surface 

*  J 

element  (do.. ). 

One  disadvantage  of  the  point  mass  model  based  on  the  ellipsoidal 
surface  as  opposed  to  a  higher  order  surface  would  be  noticed  if  geoid  un¬ 
dulations  (N)  should  be  the  predicted  quantities  in  an  adjustment  of  only 
gravity  anomalies.  It  has  been  indicated  (see  [Blaha,  1977],  pages  166-  168, 
on  which  the  results  presented  by  Needham  [1970]  are  recapitulated)  that 
large  errors  could  be  expected  in  N  if  the  locations  of  point  masses  were 
limited  only  to  one  area.  However,  much  smaller  errors  are  believed  to  oc¬ 
cur  in  the  predictions  of  geoid  undulations  or  gravity  anomalies  if  the  ad¬ 
justed  and  the  predicted  quantities  are  of  the  same  kind. 


Another  disadvantage  of  the  approach  considered  emerges  upon 
scrutinizing  the  basis  of  the  point  mass  algorithm.  If  no  initial  adjust¬ 
ment  in  terms  of  spherical  harmonics  takes  place,  the  ’’residuals"  describ¬ 
ing  the  discrepancy  between  the  observed  and  the  ellipsoidal  surfaces  will 
be  relatively  large  (they  may  reach  100m  or  more  in  some  areas).  The  role 
of  the  point  mass  adjustment  is  to  model  these  discrepancies.  However, 
since  the  point  mass  model  contains  the  spherical  approximation  and  thus 
produces  results  reliable  to  no  more  than  2  or  3  significant  digits,  this 
approach  would  be  unable  to  yield  results  even  approaching  a  decimeter  ac¬ 
curacy.  In  fact,  the  geoid  determination  could  be  affected  by  errors  of 
several  decimeters,  due  to  the  spherical  approximation  alone. 

On  the  other  hand,  a  global  adjustment  of  satellite  altimetry  in 
terms  of  fairly  high  degree  and  order  spherical  harmonic  coefficients  pro¬ 
duces  residuals  of  only  a  few  meters.  If  entered  into  a  subsequent  adjust¬ 
ment  in  terms  of  point  mass  parameters,  these  residuals  will  yield  results 
which  will  be  essentially  unaffected  by  the  spherical  approximation. 

In  summary,  upon  considering  the  effect  of  the  spherical  approxi¬ 
mation  in  the  point  mass  model,  it  becomes  clear  that  the  accuracy  require¬ 
ments  are  a  limiting  factor  in  adjustments  of  point  mass  parameters  alone. 
Indeed,  instead  of  referring  a  point  mass  adjustment  to  tie  figure  of  an  el¬ 
lipsoid,  an  introduction  of  a  higher-order  surface  as  obtained  through  a 
global  spherical  harmonic  adjustment  is  recommended.  Based  on  the  altimetry 
residuals  in  a  region  of  high-density  altimeter  data,  a  point  mass  adjust¬ 
ment  can  then  follow  in  order  to  produce  a  detailed  representation  of  the 
local  geoid. 


>.  'iM.oMlj  OXlJl.  X  lit  X  J  /All  /I  IM  A  IO<  Al  I  KAMI 
DEVELOPED  IN  SPHEROIDAL  AND  SPHERICAL  COORDINATES 


3.1  Introduction 

In  some  geodetic  applications  use  is  made  of  a  3x3  symmetric 
matrix  of  second-order  derivatives  of  the  geopotential  (W)  with  respect  to 
the  length  elements  along  local  Cartesian  coordinate  axes.  This  matrix  re¬ 
lates  the  symmetric  Marussi  tensor  to  the  five  curvature  parameters  of  the 
field  and  to  the  vertical  gradient  of  gravity  as  can  be  gathered  e.g.  from 
the  formulas  (12.162)  of  [Hotine,  1969].  One  of  the  Cartesian  axes  points 
east,  one  points  north  and  one  points  toward  the  local  zenith.  Such  an 
orientation  of  Cartesian  axes  at  a  given  point  in  space  (P)  is  linked  to 
the  equipotential  surface  W=  constant  passing  through  P.  The  directions  of 
these  local  Cartesian  axes  are  revealed  by  various  sensors  and  are  thus 
physically  meaningful. 

However,  it  is  sometimes  expedient  to  work  with  orthogonal  triads 
whose  orientation  varies  slightly  from  that  of  the  triad  just  mentioned. 

The  orientation  may  correspond  to  the  surface  U= constant  where  U  is  the 
standard  potential,  to  the  surface  a=  constant,  a  spheroid,  where  a  is  the 
third  coordinate  in  the  {u,u,a}  spheroidal  coordinate  system  described  in 
Chapter  22  of  [Hotine,  1969],  or  to  the  surface  r=  constant,  a  sphere,  where 
r  is  the  third  coordinate  in  the  {X,^,r)  spherical  coordinate  system,  X  and 
4>  being  the  geocentric  longitude  and  latitude,  respectively.  Furthermore, 
the  function  of  position  whose  second-order  derivatives  are  sought  may  be 
not  only  W  but  also  U,  T  (disturbing  potential,  T  =  W-U),  a  generalized  coor¬ 
dinate  a,  r,  etc.  Before  a  specialization  is  made,  such  a  function  will  be 
assigned  the  general  notation  F. 
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I 


1 


The  unit  vectors  along  the  local  Cartesian  axes  xx,  x2,  and  x3, 
corresponding  respectively  to  general  east,  north  and  "up"  directions,  are 
denoted  as  ^(2)  and  ^(3)  where  the  superscript  "r"  indicates  the  vec¬ 
tor's  contravariant  components  (in  Cartesian  coordinates,  covariant  and  con- 
travariant  components  are  Identical).  In  [Hotine,  1969],  abbreviated  hence¬ 
forth  as  [H],  the  above  right-handed  orthonormal  triad  is  denoted  as  Ar,  ur, 
vr,  but  for  our  purpose  it  is  advantageous  to  associate  the  i-th  unit  vector 
with  the  subscript  "(i)".  In  the  local  Cartesian  coordinates  the  triad  is 
depicted  as 

^(1)  *  f1,0,0*  **•  Poising  in  general  east, 

^(2)  =  •**  Pointing  in  general  north, 

^(3)  =  f0’0*1)  •••  pointing  in  general  "up". 

The  desired  second-order  derivatives  of  F  with  respect  to  the 
length  elements  along  the  local  Cartesian  axes  xl,  x2,  x3  can  themselves  be 
viewed  as  scalar  functions  of  position;  any  of  them  can  be  expressed  as 

32F/3x13xj  =  (32F/3xr3xS)£[ij^J.j,  (3.1) 

where  use  is  made  of  the  summation  convention  for  repeating  indices  and 
where  the  symmetric  matrix  composed  of  32F/3x  3x  for  r,s= 1,2,3  is 


{32F/3xr3xs} 


32F/3x13x1 

32F/3xx3x2 

32F/3xx3x3 


32F/3x*3x2 
32F/3x23x2 
32F  3x23x3 


32F/3xx3x3 

32F/3x23x3 

32F/3x*3x3 
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As  in  this  example,  {  }  will  indicate  a  matrix  where  either  of  the  indices 
assumes  the  values  1,2,3.  Another  notation  convention  which  will  be  adhered 
to  is  the  replacement  of  dx1  by  ds^j  as  a  length  element  along  the  axis  x^. 
Since  in  Cartesian  coordinates  32F/3xr3xs  are  identical  to  the  second  co¬ 
variant  derivatives  Fr$,  (3.1)  can  be  rewritten  as 

32F/3s(i)3s(j)  =  Frs*(i)£(j>'  (3*2) 

But  this  is  a  tensor  invariant  which,  as  its  name  indicates,  has  the  same 
value  regardless  of  the  coordinate  system  in  which  the  given  vectors 
and  the  covariant  derivatives  Fr$  are  expressed.  The  six  distinct  values 
obtained  from  (3.2)  with  (i),(j)  =  1,2,3  correspond  to  the  left-hand  sides 
of  the  six  equations  in  (12.162)  of  [h],  provided  our  notations  F, 

£(2).  £(3)  a^e  made  to  coincide  with  the  notations  N,  X,  p,  v  of  [H],  re¬ 
spectively. 

Later  on  when  the  final  results  are  listed,  the  more  usual  geo¬ 
detic  conventions  will  be  adopted  where  the  first  axis  (x)  corresponds  to 
the  northerly  direction,  the  second  axis  (y)  corresponds  to  the  easterly 
direction  and  the  third  axis  (z)  corresponds  to  the 
the  following  transliterations  will  take  place: 

32F/3s(l)3S(i)  =  32F/3y2  =  Fyy, 

32F/3s^j3S^j  =  32F/3y3x  =  32F/3x3y  5 

32F/3S(1)3S(3)  =  32F/3y3z  =  Fyz> 


'up"  direction.  Thus, 


Fxy* 


(3.3) 
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32F73S(2)3s(2)  =  32F/3x2  3  F^, 

32F/3S(2)3S(3)  3  32F/3x3z  3  Fxj,, 

32F/3S(3)3S(3)  =  32F/3z2  3  F^. 

As  indicated  at  the  outset,  more  than  one  orthonormal  triad  ema¬ 
nating  from  P  will  be  considered.  In  fact,  four  right-handed  orthonormal 
triads  will  be  utilized  where  the  index  in  parentheses  identifies  not  only 
the  particular  unit  vector  of  a  triad  but  also  the  triad  itself.  Upon  de¬ 
noting  the  geocenter  by  0,  the  following  description  applies  with  regard 
to  the  third  vector  of  each  triad: 

£^3j  ...  vector  orthogonal  to  a  sphere  centered  at  0, 
vector  orthogonal  to  a  spheroid  centered  at  0, 

iLj...  vector  orthogonal  to  an  equipotential  surface 
1  '  (spherop)  in  the  standard  gravity  field,  cen¬ 

tered  at  0, 


£(3\  ...  a  general  vector  coplanar  with  the  other  three 
'  '  vectors  above. 

This  situation  is  illustrated  in  Figure  3.1  with  some  additional  notations; 
for  example,  $  indicates  the  direction  of  the  line  of  force  at  P  in  the 
standard  gravity  field  (in  agreement  with  [H],  the  word  "standard"  is  used 
in  lieu  of  "normal").  From  the  figure  one  reads 

♦  =  4>  +  »  (3.4a) 

$  =  $  +  6Z  *  $  +  Sl  +  52;  (3.4b) 


-15- 


sometimes  the  notation 


6  =  &x  +  62  (3.4c) 

will  be  utilized.  Figure  3.1  depicts  also  whl1e  tfj)  ~  ^(l)  *  ^(i')  = 

go  into  the  plane  of  the  paper  at  P. 


Schematic  representation  of  four  kinds  of  n'ght-handed 
orthonormal  triads,  i,  i',  i"  and  i 


The  triad  i  is  obtained  from  the  triad  i'  by  means  of  the  rotation 
Rj(-t),  giving  rise  to  the  following  vector  equations  valid  in  any  coor¬ 
dinates: 


tT 

1')  ’ 

cost  £^2' )  -  sinx  j 
sinx  ^2')  +  C0ST  ^(3 ' ) 


j 


(3.5) 
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in  the  matrix  notations,  these  equations  are  represented  by 


U(T)}  =  Rl<-T)  ll(V)}-  (3.5*) 

The  special  cases  of  interest  are  the  following: 

t  =  §2  •••  triad  i  becomes  triad  i"  , 

x  =  0  ...  triad  i  becomes  triad  i‘  ,  >  (3.6) 

x  *  ...  triad  i  becomes  triad  i  .  J 

Starting  with  equation  £-2),  we  keep  in  mind  that  tensor  or  vec¬ 
tor  equations  are  valid  in  any  coordinates,  such  as  local  Cartesian,  earth- 
fixed  Cartesian  (which  need  not  be  mentioned  again),  spheroidal,  spherical, 
etc.  Only  two  kinds  of  coordinate  systems  will  be  utilized  here  in  order 
to  develop  expressions  such  as  (3.2):  the  spheroidal  coordinate  system  and 
the  spherical  coordinate  system.  The  former  is  identified  by  primes;  we 
thus  consider  {x'1,  x'2,  x'3}  =  {u.u.a},  Vr,  l},,  FJ,  ,  etc.  In  the  latter 
the  primes  are  absent  and  the  corresponding  notations  are  (x1,  x2,  x3}  = 
{A,4>,r},  lr,  lr,  Frs,  etc.  That  all  the  expressions  so  far  have  been  writ¬ 
ten  without  primes  does  not  mean  that  they  will  be  considered  only  in  sphe¬ 
rical  coordinates.  For  example,  (3.2)  can  equally  well  be  written  as 

8!F'/as(i)3s(j)  .  ;  (3.7) 

since  F'=  F,  the  prime  on  the  left-hand  side  serves  only  to  indicate  in 
which  coordinate  system  the  computation  Is  being  made.  Similarly,  (3.5) 
can  be  written  in  spheroidal  coordinates,  i.e.,  every  lr  can  be  replaced 
by  Vr. 
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The  desired  second-order  derivatives  of  F  will  first  be  developed 
in  spheroidal  coordinates.  The  form  of  the  results  will  correspond  to 
(3.7)  in  which,  however,  use  will  be  made  of  a  more  general  triad  i,  i.e., 

92F73S(T)8S(J)  =  F^(T)£(f)  .  (3.8) 

The  other  cases,  including  (3.7),  can  be  deduced  from  this  result  upon 
specializing  x  as  in  (3.6). 

Subsequently,  the  corresponding  expressions  will  be  obtained  in 
spherical  coordinates.  However,  this  will  not  be  done  directly  using  the 
matrix  tensor,  the  associated  metric  tensor  and  the  Christoffel  symbols  in 
these  coordinates,  but  indirectly  using  the  already  known  Christoffel  sym¬ 
bols  in  spheroidal  coordinates  and  several  differential  relations  obtained 
from  the  geometric  properties  of  the  two  systems.  On  one  hand,  this  could 
be  considered  merely  an  academic  exercise  since  the  results  can  be  obtained 
much  faster  in  spherical  coordinates  by  direct  means.  On  the  other  hand, 
such  an  approach  could  prove  useful  if  one  wanted  to  proceed  from  the  sphe¬ 
roidal  (or  some  other)  coordinate  system  to  a  more  complicated  system  with¬ 
out  utilizing  the  metric,  the  Christoffel  symbols,  etc.,  in  this  new  system. 
(Although  the  associated  metric  tensor  in  this  new  system  can  be  derived 
through  the  same  differential  relations  as  mentioned  above  and  although  the 
metric  tensor  as  well  as  the  Christoffel  symbols  can  be  derived  subsequently, 
such  a  process  might  not  be  always  shorter  than  the  one  discussed.)  Be  that 
as  it  may,  the  link  between  the  two  approaches  can  also  serve  as  an  inde¬ 
pendent  verification  of  the  results  obtained  separately  in  one  and  the 
other  coordinate  system. 


-18- 


In  the  next  step,  the  development  will  be  carried  out  directly  in 
the  spherical  coordinate  system,  leading  to  the  results  of  type  (3.2)  but 
with  a  more  general  triad  i,  i.e., 

32F/3s(T)3s(j}  =  Frs£(T)£(T)-  (3‘9) 

The  specializations  indicated  in  (3.6)  will  again  encompass  the  other  cases 
considered,  including  (3.2).  The  approach  of  [H]  adopted  in  this  study  will 
be  seen  to  lead  to  the  desired  results  much  faster  than  more  conventional 
approaches  in  which  advantage  is  not  taken  of  tensors  and  their  properties. 

In  the  closing  part  of  this  chapter,  the  practical  applicability 
of  the  results  obtained  in  spherical  coordinates  will  be  discussed.  As  an 
example,  three  of  the  six  distinct  second-order  derivatives  of  W  can  rou¬ 
tinely  be  measured  as  one  vertical  and  two  horizontal  gradients  of  gravity. 
Upon  using  the  derived  results,  formulas  for  the  second-order  derivatives 
of  T  will  be  established  making  it  possible  to  model,  among  other  things, 
the  "anomalous  gradients  of  gravity".  Advantage  will  be  taken  of  the  La¬ 
place  equation  AT  =  0  which  will  eliminate  the  need  for  the  second  derivatives 
of  the  associated  Legendre  functions  with  respect  to  3>,  should  these  be  a 
part  of  the  explicit  model.  The  desired  second-order  derivatives  of  T  will 
be  expressed  for  a  general  point  in  space,  whether  on  the  earth's  surface 
or  above  it  (e.g.  at  aircraft  or  satellite  altitudes).  The  results  will  be 
further  generalized  by  being  made  applicable  to  a  rotating  field  as  well 
as  a  nonrotating  field. 
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3.2  Derivations  in  Spheroidal  Coordinates 

In  order  to  evaluate  various  quantities  in  spheroidal  coordinates 
{to,u,a}  the  values  of  these  coordinates  should  be  known.  If,  for  example, 
the  earth-fixed  coordinates  are  given,  the  spherical  coordinates  {A,<f>,r} 
can  be  computed  immediately  from 

tgX  =  Y/X,  I 

tg*  =  Z/(X2  +  Y2)^  >  >  (3.10) 

r  =  (X2  +  Y2  +  Z2)%  .  J 

When  working  with  spheroidal  coordinates  the  situation  is  more 
complicated.  The  bulk  of  the  pertinent  geometry  is  presented  on  pages 
187  -  190  of  [H].  We  shall  now  review  some  of  the  basic  geometric  relation¬ 
ships  with  the  aid  of  Figure  3.2.  The  absolute  space  constant  (ae),  con¬ 
sidered  known,  gives  rise  to  a  family  of  spheroids  (rotational  ellipsoids) 

one  of  which  passes  through  P.  The  eccentricity  (e),  the  semi-major  axis 
(a)  and  the  semi-minor  axis  (b)  of  the  spheroid  a= constant  passing  through 
P  are  related  as  follows: 

e  =  sina  ,  (3.11a) 

b  =  aO-e2)*5  =a  cosa  .  (3.11b) 

Except  for  the  (well-known)  last  equality  in  (3.12a)  below,  all  the  relations 
in  (3.12)  can  be  read  directly  from  Figure  3.2: 
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PR  =  r  sin<J>  =  b  sin  u  =  v  cos2a  sin<{>  , 


(3.12a) 


OR  =  r  cos<J>  =  a  cos  u  =  v  cos<|>  ;  (3.12b) 

It  thus  follows  that 

tg4>  *  cosa  tg  u  =  cos2a  tg4»  .  (3.13) 

From  the  same  figure  and  from  the  definition  of  an  ellipse  we  have 
d|  =  [r  cos?  +  (ae)]2  +  (r  sin?)2  , 

d|  =  [r  cos?  -  (ae)]2  +  (r  sin?)2  , 
a  =  ^(dj  +  dg). 


The  spheroidal  coordinates  are  thus  computed  from 


to  *  X  , 
sina  =  (ae)/a, 
tgu  =  tg<|>/cosa  , 


(3.14) 


where,  in  addition  to  (ae),  also  X  and  <J>  are  known  beforehand  (see  e.g. 
equations  3.10).  Upon  consulting  (3.13)  and  (3.4a)  or  Figure  3.2,  one 
can  write 


tg<j>  *  tg  u/cosa  *  tg<f>/cos2a  ,  (3.15) 

\  (3.16) 
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F,0  =  0Fo  =  (ae)  =  (a  e  )  =  absolute  constant 
1  2  oo 

AjO  =  0A2  =  OP'  -  a 


A  few  elements  of  the  geometry  pertinent 
to  the  spheroidal  coordinate  system 


From  (22.25)  of  [h],  the  metric  in  spheroidal  coordinates  is 
ds2  =  (a2cos2u)du>2  +  (v2cos2a)du2  +  (v2cotg2a)da2  . 

This  results  in  the  following  metric  and  associated  metric  tensors: 

{gj,$}  =  diag.  {a2cos2u,  v2cos2a,  v2cotg2a)  ,  (3.17a) 

{g'rs}  =  diag.  (l/(a2cos2u) ,  l/(v2cos2a),  l/(v2cotg2a)>.  (3.17b) 

In  a  triply  orthogonal  coordinate  system  we  have  along  the  coordinate  lines: 
dS(i)  =  (g^J^dx1  ,  etc.  When  introduced  in  the  general  formula  £r  =  dxr/ds, 
this  yields  =  {(g11)*5.  0,0},  etc.  Applying  this  outcome  to  the 


(triply  orthogonal)  spheroidal  coordinate  system  and  adhering  to  the 
earlier  conventions,  we  deduce 


■fcjj.)  =  (l/(a  cos  u),  0  ,  0  }  , 

2')  =  {0  •  1/ ( vcosot)  ,  0  )  ,  > 
•d(3')  =  (0,0,  -l/(vcotga)>  .  ^ 


(3.18) 


The  minus  sign  in  the  last  relation  above  is  due  to  a  increasing  inward; 
by  contrast,  the  third  vector  of  our  orthonormal  triad  is  an  outward-drawn 
normal.  To  reconcile  this  difference  we  have  adopted  -da/ds^ij  as  the 
vector's  third  component.  Taking  (3.12b)  into  consideration,  in  agreement 
with  (3.5)  we  write 


U/(vcos<J>)  ,0,0), 

(0  ,  cosx/(vcosa)  ,  sinx/(vcotga))  , 

(0  ,  sinx/(vcosa)  ,-cosx/(vcotga)}  .  > 


(3.19) 


The  Christoff el  symbols  of  the  second  kind,  needed  in  forming 
second  covariant  derivatives  of  F,  are  defined  in  any  coordinates  as 


r£s  =  9tktrs»  *3  •  (3.20a) 

[rs,  k]  =  »s(3grk/3xs  +  3gsk/axr  -  ag^/ax1).  (3.20b) 
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The  symbols  r*s  are  symmetric  in  the  lower  indices  and  as  such,  are  often 
denoted  by  the  symbol  {£  }  in  tensor  calculus.  However,  in  the  present 
context  the  symbolism  of  [h]  will  be  utilized.  The  above  formulas  are  ap¬ 
plied  to  our  system  in  a  straightforward  fashion.  In  the  process,  the  fol¬ 
lowing  differential  relations  are  utilized: 

3a/ 3w  =  0,  3a/3u  =  0,  3a/3a  =  -a  cotga  , 

3v/3w  =  0,  3v/3u  =  (a2/v)tg2a  sin  u  cos  u,  3v/3a  =  v(-cotga+tgasin2<|>). 


After  some  manipulation  we  obtain  essentially  (22.38)  and  (22.39)  of  [H] 
which,  in  this  reference,  were  derived  somewhat  differently.  The  nonzero 
Christoffel  symbols  (symmetric  in  the  lower  indices)  are  listed  as 


=  -tg  u  =  -cosa  tg<J>  , 

ri3  =  "cot9a  » 

,  2 

rxl  =  sin<p  cost))/cosa  , 

r22  =  S1'na  *9®  sin<^  cos<*)  * 

-  -a2/(v2sina  cosa)  , 


(3.21) 


-sin«f>  cos^/cosa  ; 


rjJ  =  tga  cos24>  , 

*22  =  (a2/v2)tga  , 

=  sina  tga  sin<j>  cos$  , 

3 

I33  =  -cotga  -  a2/(v2  sina  cosa)  . 

The  general  formula  giving  the  second  covariant  derivatives 
(symmetric)  of  F  reads 

Frs  =  "  r£sFt  »  Ft  =  3F/9X1  .  (3.22) 

In  spheroidal  coordinates  this  results  in 

Fji  =  32F/3w2-  ( si n<p  cost})/ cosa) 3F/ 3u  -  tga  cos2<}>  3F/3a  , 

Fj2  *  32F/3(i)3u  +  cosa  tg<p  3F/3u  ,  j 

Fjs  =  32F/3(o3a+ cotga  3F/3u  , 

F22  =  92f/9ij2  -  sina  tga  sin<j>  cos4>  3F/3u  -  (a2/v2)tga  3F/3a  , 

F23  =  92F/9u9a+  [a2/(v2sina  cosa)]  3F/3u  -  sina  tga  sin<p  cos<i>  3F/3a  , 

Fj3  =  32F/3a2  +  (sin<J»  cos<(>/cosa)3F/3u  +  [cotga+  a2/(v2sina  cosa)]  3F/3a  . 

We  can  proceed  according  to  the  formula  (3.8),  symmetric  in  (T), 
(j),  in  which  the  values  (3.19)  and  (3.23)  are  utilized.  The  results  are 


)  (3.23) 
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3 2F */3s^Y) 3s (Y)  =  (l/v2cos24>)92F/9u)2  -  [tg<j>/(v2cosa)]  9F/9u 

-  [l/(v2cotgoi)]  9F/9a, 

92F73SY)9s(2)  =  [cosx/(v2cosacos<}>)]  32F/9u>9u 
+  [sinx/(v2cotgacos<j>)]  92F/9a)9a 
+  [sin(<j>+x)/(v2cos2<J))]  9F/9u), 

92F'/3s =  [sinr/(v2cosacos<j>)]  92F/9oj9u 

-  [cosx/(v2cotgacos<J>)]  92F/9w9a 

-  [cos(<ji+x)/(v2cos2<j>)]  9F/9tu, 

92F,/3s^j9s^  =  [cos2t/(v2COS2a)]  92F/9u2 

+  [sin  2x/(v2COSacotga)]  92F/9u9a 
+  [sin2x/(v2cotg2a)]  92F/9a2 

-  [l/(v2cosacotg2a)]  [cos  2x  «  sinc|>cos(|> 

-  sin  2x  .  a2/(v2sin2a)] 9F/9u 

-  [l/(v2cotg3a)]  [cos  2x  .  a2/(v2sin2a) 

+  sin  2x  «  sin4>cos4>  -  sin2xcotg2a]9F/9a, 

92F'/3s^2j3s^yj  =  ^[sin  2x/(v2cos2a)] 92F/9u2 

-  [cos  2x/{v2cosacotga)]  92F/9u9a 

-  *s[sin  2x/(v2cotg2a)]  92F/9a2 

-  [l/(v2cosacotg2a)]  [cos  2x «  a2/(v2sin2a) 

+  sin  2x»  sin<|)COS4)]9F/9u 

+  [l/(v2cotg3a)]  [cos  2x «  sin<f>cos<f> 

-  sin  2x«  a2/(v2sin2a) 

-  h  sin  2x  .  cotg2a]9F/9a, 


(3.24) 


32F78S(j)3s^  =  [sin2T/(v2cos2ot)j  32F/3u2 

-  [sin  2t/(v2cosacotga)]32F/3u3a 
+  [cos2t/(v2cotg2a)]  32F/3a2 

+  [l/(v2cosacotg2a)]  [cos  2i  •  sin<j>cos<f> 

-  sin  2t «  a2/(v2sin2a)] 3F/3u 

+  [l/(v2cotg3a)] [cos  2x «  a2/(v2sin2a) 

+  cos2t  cotg2ct+sin  2t  »  sin<j>cos<{»]3F/3a. 


By  setting  r  to  »  0  »  or  -6^  we  obtain  the  desired  second- 
order  derivatives  of  F  for  the  triads  1"  ,  i1  ,  or  i  ,  respectively,  as 
indicated  in  (3.6).  The  choice  r=0  results  in  particularly  simple  formulas 
they  can  be,  of  course,  obtained  also  more  directly  upon  replacing  in  (3.8) 
the  triad  (3.19)  by  the  triad  (3.18).  The  expression  for  the  Laplacian  of 
F  in  spheroidal  coordinates  (AF' )  is  obtained  as  the  trace  of  the  symmetric 
matrix  composed  of  the  elements  in  (3.24)  regardless  of  the  choice  of  t, 
i.e. , 

AF'  =  [l/(v2C0S2<j>)]32F/3w2  +  [l/(v2cos2a)]  32F/3u2  +  [l/(v2cotg2a)J  32F/3a2 

-  [tg«J>/  (v2co$a)]  3F/3u ; 

a  similar  form  appeared  in  (22.41)  of  [h].  Since  the  Laplacian  is  a  tensor 
invariant  the  same  value  must  be  obtained  in  any  coordinates  (spheroidal, 
spherical,  Cartesian,  etc.). 
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The  computation  of  6^  is  an  easy  matter  as  was  already  shown  in 
(3.16).  On  the  other  hand,  the  computation  of  62  involves  a  formula  ex¬ 
pressing  U  (standard  potential,  in  Chapter  23  of  [H]  denoted  as  -W)  at  P. 
The  main  task  consists  in  projecting  the  vector  IV  (=3U/3x'r)  onto  the 
triad  (3.18),  the  projections  being  respectively 


3U7as(1>)  =  U7jf()  =  0  , 

3U7as(2,j  =  ur^(2. )  =  [l/(vcoso)]au/3u  , 
3U73s(3,j  =  li'Uy)  =  -[l/(vcotga)]3U/3a  . 


(3.25) 


The  expression  for  U  is  given  in  (23.13)  of  [H]  essentially  as 

(3.26) 

U  -  GMa/(ae)  +  GA2QQ2(i  cotga)  P2(sin  u)  +  [<L2a2/3  -  aj2a2P2(sin  u)/3], 

where  GM  is  the  product  of  the  gravitational  constant  and  the  earth's  mass, 
(5  is  the  angular  velocity  of  the  earth's  rotation,  P2(sinu)  is  the  Legendre 
function  in  argument  sin  u  and,  according  to  (22.52)  and  (23.12)  of  [H], 


Q2(i  cotga)  =  W  (a+  3a  cotg2a  -  3cotga)  ,  (3.27a) 

iGA20  =  ~(2/3)u2a2/ [3cotga0  -  a0(l  + 3c0tg2a0)]  ,  (3.27b) 

where  aQ  ,  aQ  correspond  to  a  given  base  equipotential  spheroid.  Outside 
this  particular  spheroid  the  equipotential  surfaces  U= constant  are  not 
spheroids.  It  is  to  be  noted  that  if  the  standard  field  should  be  nonrota¬ 
ting,  such  as  considered  in  some  satellite  applications,  the  terms  in 
(3.26)  containing  S32  explicitly  would  disappear  (the  numerical  value  of 
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the  constant  iGA^  in3.27bwould  remain  unchanged).  In  this  case,  none  of 
the  modified  surfaces  U  *  constant  would  ever  coincide  with  a  spheroid. 

The  angle  62  would  be  computed  exactly  as  in  the  present  approach  except 
that  the  explicit  symbols  w2  in  the  pertinent  expressions  would  be  replaced 
by  zeros. 

Carr  ing  out  the  indicated  operations,  from  (3.25)  we  obtain 
3U73S(2,)  =  03  sin  u  cosu/  (vcosa)]  (GA20Q2(i  cot9a)  "  s2a2/3])  ;  (3.28a) 

on  the  base  spheroid  the  expression  in  the  second  brackets  becomes  zero  (in 
the  rotating  field).  Further  we  have 

3U73S(3,j  =  -[l/(vcotga)]  [J  + LP2(sinu)]  ,  (3.28b) 

where,  similar  to  (23.31)  of  [H], 

J  *  GM/(ae)  -  (2/3)G2a2cotga  , 

L  s  -iGA2g(-2  -  3cotg2a  +  3a  cotga+3a  cotg3a) 

+ (2/3)u>2a2cotga  . 

Since  both  (3.28a),  (3.28b)  are  in  general  negative  outside  the  equipoten- 
tial  spheroid  (rotating  field),  the  desired  angle  <$2  1S  obtained  from 

tgfi2  =  (-3U73s(2,j)/(-3U73Sp,))  ,  (3.29) 

and  it  is  in  general  positive  (at  the  poles  or  on  the  equator  it  is  zero). 
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In  the  nonrotating  field  the  value  of  6^  would  be  usually  negative;  it 
can  be  shown  that  this  sign  would  start  changing  at  high  altitudes,  ap¬ 
proximately  at  a^/3  (about  2.1xl06m)  above  the  base  spheroid. 

Denoting  the  standard  gravity  at  P  by  y  »  its  normal  component 
(positive  inward  along  -ds^.j)  by  Yn  and  its  meridian  component  (positive 
north  along  ds^.^)  byy^,  we  have 


Yn  =  -SU'/Ss^i^  ♦  (3.30a) 

Ym  =  3U73s(2.j  ,  (3.30b) 

Y  »  (Yj  +  V*)*5  ,  (3.30c) 

and  from  (3.29), 

tg  62  =  -Ym/Yn  •  (3.31) 


We  notice  that  the  Somigliana  formula  appearing  in  (23.24)  of  [h]  gives  Yn 
at  any  point  in  space  (on  the  base  spheroid  it  follows  from  3. 28a  that  Ym=0 
and  Y«Yn).  However,  Yn  given  by  that  formula  is  computed  from  Yg  and  Yp 
representing  the  standard  gravity  on  the  equator  and  at  the  poles  of  the 
coordinate  spheroid  passing  through  the  point  in  question  (if  it  is  not 
the  base  spheroid  such  a  spheroid  is  not  an  equipotential  surface).  In 
the  applications  involving  the  base  spheroid  the  rigorous  Somigliana  for¬ 
mula  is  quite  suitable.  However,  since  in  general  for  each  point  in  space 
the  corresponding  Ye  and  Yp  would  have  to  be  computed  first,  the  present 
approach  giving  Yp  in  (3.30a)  in  conjunction  with  (3.28b)  appears  to  be 
more  practical.  This  assertion  would  hold  all  the  more  if  the  nonrotating 
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field  were  considered  instead  (  such  a  case  all  the  formulas  of  the 
present  approach  would  be  simplified  upon  replacing  the  explicit  symbols 
u>2  by  zeros). 

Having  obtained  62  »  the  desired  derivatives  can  further  be  re¬ 
fined  so  as  to  refer  to  an  i  system  which  characterizes  the  actual  gra¬ 
vity  field.  This  system,  which  could  be  termed  as  a  system  of  the  local 
vertical,  differs  from  the  i"  system  by  two  rotations,  each  corresponding 
to  one  component  of  the  deflection  of  the  vertical  at  P.  In  general,  the 
second  covariant  derivatives  in  any  two  coordinate  systems  are  connected 
by 

{F.j}  =  {ax,,r/3xi}T{F^sH3x"s/8x'j}  ; 

the  indices  are  unimportant  here  since  we  are  dealing  with  matrix  expres¬ 
sions.  If  the  two  systems  are  Cartesian  —  and  thus  the  Christoffel  sym¬ 
bols  vanish  —  one  can  write 

{32F/3x13xJ}  =  CT{32F/3x,,r3x,,S}C  ,  (3.32) 

CT  =  {3x1/3x"r}  ,  {dX1 }  =  CT{dx"r}  . 

In  our  case  CT  represents  the  following  rotation  matrix: 


(3.33) 
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where  £  ,  n  are  the  north  and  east  components  of  the  deflection  of  the 
vertical  at  P.  If  (3.33)  is  utilized  in  (3.32)  and  the  previous  conven¬ 
tions  are  adhered  to  (dx1  written  as  ds^  ,  F  written  as  F*)»  it  follows 
that 


32F,/3s^!)a$(i) 

32F'/3s^2^3s^2j 

32F'/3S(i)3s(3) 

32F'/3s^2j9s^2j 

32F73S(2)3S(3) 

32F73S(3)3S(3) 


-  3 2F 73s ^ 3s ^ j  -  2r)  32F79s^n^3s^3u^  , 

-  32F  73s  ^  3s^  2"  j 

-  £  3  F  /3s^jhj  Ss^n )  -  n  3  F  /3s^2h  )3s(3"  )  ’ 

~  32F 73s jjn ^ 3s ^  +  £  32F73s^,i ^3s^2» j 
n[92F73s^3ii  j3s^3„  j  -  32F73s^,i  j3s^„  j]  , 

-  32F  /3s ^ 2»  j ^ 2" )  ~  32F  /3s^2"  )^s(3" )  ’ 

-  3  F  /  8s  ^  2*« )  ( 3"  )  ~  ?[32F  /^s(3")^s(  3") 

-  3  F  /3s ^ 2" ) ^S ^ 2" ) J ^  h  3  F  /3s^jh j3s^2»j  » 

-  32F73S^3„j3s^3„^  +  2^  32F 73S| 2m ^s^ii ^ 

+  2n  32F 73s^m jSs^,, ^  . 


(3.34) 


These  relationships  will  be  used  also  later  when  dealing  with  spherical 
coordinates,  in  which  case  the  primes  associated  with  F  will  be  dropped. 
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3.3  Transformation  of  the  Results  from 
Spheroidal  to  Spherical  Coordinates 


The  expressions  in  (3.8)  and  (3.9)  being  tensor  invariants  yield 
the  same  values  in  any  coordinates,  i.e.. 


32F/3s(i)3s(j)  =  Frs 


fr-  fs- 

(i)(j) 


f;s 


pijc.  p  i  s. 

€dr(j) 


(3.35) 


In  spheroidal  coordinates  the  last  relation  resulted  in  (3.24)  featuring 
the  partial  derivatives  32F/3w2  ,  ...  ,  3F/3a  .  It  is  now  required  that 
the  second-order  derivatives  (3.35)  be  written  in  terms  of  32F/3A2  ....  , 
3F/3r  .  However,  this  task  is  to  be  accomplished  without  resorting  to  a 
development  in  spherical  coordinates  which  is  the  subject  of  the  next  sec¬ 
tion.  This  amounts  to  avoiding  the  formulation  of  the  Christoffel  symbols 
(rj  )  in  spherical  coordinates  --  or  in  some  other  coordinates  if  we  wish 
to  generalize  this  approach  —  and  taking  instead  advantage  of  r'*s  al¬ 
ready  found.  The  price  to  pay  is  the  need  for  certain  differential  rela¬ 
tionships  between  the  two  systems. 

The  straightforward  approach  adopted  in  this  derivation  proceeds 

from 

F;s  =  32F/ax,r3x,S  -  r'JsF'  ,  (3.36) 

where 

FJ.  =  (3xk/3x,t:)Fk  ;  Fj.  =  3F/3x't,  Fk  n  3F/3xk  .  (3.37) 
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The  first  member  on  the  right-hand  side  of  (3.36)  is  developed  as 

92F/9x,r9x‘S  =  9FJ./9x,S  =  9  [(9xk/9x,r)Fk]/9x'S 
=  (9xP/9x'r)(9xq/9x,S)92F/9xP9xq 

+  Fk9(9xk/9x,r)/9x'S  , 

where  the  presence  of  the  second  term  stems  from  the  fact  that  in  general 
coordinates  a  second  derivative  is  not  a  tensor.  If  this  expression  as 
well  as  (3.37)  are  substituted  in  (3.36)  and  if  the  new  equation  is  con¬ 
tracted  with  as  indicated  in  (3.35),  one  obtains 

(i)  (j) 

92F/9s^j9s^jj  =  (92F/9xr9xS)£^f^  +  [9(9xk/9x,r)/9x,S 

-  r^s9xk/9x,t]f^)f[^)9F/9xk  ,  (3.38) 

where  use  has  been  made  of  the  general  transformation 

lm  =  (9xm/9x,n)rn  .  (3.39) 

Equation  (3.38)  is  the  desired  formula  in  which  all  the  derivatives  involve 
only  the  spherical  coordinates,  yet  the  Christoffel  symbols  in  these  coor¬ 
dinates  are  absent. 

As  a  matter  of  interest  we  may  equate  (3.38)  with  the  middle 
member  of  (3.35)  developed  according  to  (3.22).  It  then  follows  that  the 
second  term  on  the  right-hand  side  of  (3.38)  must  be  equal  to 
"rrs^(T)^(j’)F|<;  for  a"y  •£(-{) »  and  F.  Thus  it  must  also  hold  that 
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rpq  =  Cr'Js3xk/3x,t  "  3(9xk/3x'r)/3x,s](5x,r/Sxp)Ox's/9xq). 


(3.40) 


The  presence  of  the  second  term  within  the  brackets  in  (3.40)  is  due  to  the 
Christoffel  symbols  not  being  tensors.  When  the  Christoffel  symbols  on  the 
left-hand  side  of  (3.40)  were  computed  in  this  manner,  the  agreement  with 
those  derived  directly  in  spherical  coordinates  was  perfect.  This  confirms 
the  formula  (3.38)  of  which  (3.40)  is  a  by-product. 

An  important  step  in  the  present  development  is  the  formation  of 
the  matrix  {3xk/3x'r}  which,  in  (3.38),  serves  to  obtain  '  ,  t*  from 

CO  (7) 

their  counterparts  in  spheroidal  coordinates  as  is  shown  in  (3.39).  Fur¬ 
thermore,  the  elements  of  this  matrix  have  to  be  differentiated  with  re¬ 
spect  to  x,s.  In  addition  to  the  formulas  (3.11)  -  (3.16),  some  of  the 


more  important  relationships  needed  in  these  tasks  are: 

sin61  =  siVasin^coslF  =  (v/r)sin2asin<j>cos<t>  ,  (3.41a) 
cosfij  =  a2/(vr)  ,  (3.41b) 
sin  u  =  (v/a)cosasin<|>  ,  (3.41c) 
cos  u  =  (v/a)cos<J>  ,  (3.41d) 
v  =  a(l  -  sin2acos2u)Vcosa  =  a/(l  -  sin2asin24>) 2  ,  (3.41e) 


1  -  sin2asin2<J>  =  cos2a  +  sin2acos2<J>  =  a2/v2  =  (r/vJcosSj  .  (3.41f) 
After  some  manipulation  we  deduce 
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__  .* 


(3.42) 


9A/3u)  9A/9U  3A/9a 

— 

1  0 

— 

0 

3(j)/3(i)  d<p/du  94>/9a 

= 

0  (v/r)cosacos<5^ 

-(v/r)cotgasin6j 

9r/9ui  9r/9u  9r/9a 

0  -vcosasinej 

-vcotgacosfij 

Since  this  matrix  can  be  decomposed  into  a  product  of  three  matrices  of 
which  two  are  diagonal  and  one  is  a  rotation  matrix,  (Sx'^/Dx^)  can  be 
formed  without  any  matrix  inversions.  However,  the  latter  transformation 
matrix  is  not  needed  in  this  presentation. 

By  applying  (3.39)  and  (3.42),  the  vectors  appearing  in  (3.19) 
can  now  be  expressed  in  spherical  coordinates  as 


=  ( 1/ ( r  cos<p)  ,0,0}, 


■er_  =  {0  ,  cos(<S1+x)/r  ,  -sin(6,+T)}  ,  V 
(2)  1  1  ( 


=  (0  ,  sin(6,+x)/r  ,  cos(6.+t)}  . 

(3)  1  1 


(3.43) 


In  order  to  differentiate  the  elements  of  (3.42)  with  respect  to 
x,s  the  following  auxiliary  relations  are  derived: 


9<j>/  9u) =  0, 
9v/9<d  =  0, 
3<S  j  /  9cu  =  0, 


9<f>/9u  =  r  cos61/(vcos«) ,  9<{>/9a=r  sin6^/(vsinacosa) ; 

9v/9u  =  r  sinfij/cosa,  9v/9a  =  vtga-  r  cos6j/(sinacosa) ; 
96j/9u  =  r  cos6j/(vcosa)  -  (v/r)cosacos61  , 

96j/3a=r  sin6j/(vsinacosa)  +  (v/r)cotgasin6^  . 


After  some  manipulation,  3(3x  /3x'  )/3x 1  is  expressed  in  the  following 
(symmetric)  pairs  "(r,s)": 


(1.1)  ...  3(1,0, 0}/3w  =  {0,0,0}  , 

(1.2)  ...  {0,0,0}  , 

(1,3)  ...  {0,0,0}  , 

(2.2)  ...  {0,  (v2cos2a/r2)sin  26j  ,  -r +  (v2cos2a/r)cos261}  , 

(2.3)  ...  {0,  -l/sina+  (v2/r2)cosacotgacos  26^  , 

-{v2/r)cosacotgasin6jCOs6j}, 

(3.3)  ...  {0,  (v/rjcotg^sinfij  -  (v2/r2)cotg2asin  2^  , 

r/sin2a+vcotg2acos6j  +  (v2/r)cotg2asin26j}. 


t  k  t 

The  next  step  involves  the  formation  of  I”rs3x  /3x'  ,  all  the 
elements  being  already  known.  The  straightforward  algebra  produces  the 
following  results  listed  again  according  to  the  pairs  "(r,s)": 


(1.1)  ...  {0,  sin<fcos$,  -r  cos2$}  , 

(1.2)  ...  {-cosatg<f>,  0,  0}  , 

(1.3)  . . .  {-cotg<*»  0,  0}  , 

(2.2)  ...  {0,  0,  -r}  , 

(2.3)  ...  {0,  -1/sina,  0}  , 


(3,3)  ...  {0,  (v/rjcotg^sindj. 


r/sin2a  +  vcotg2acos61 }  . 


This,  together  with  the  preceding  results,  leads  to  the  expressions  for 
[3(3xk/3x'r)/3x's  -  r,t  3xk/3x't]  listed  below  again  according  to  the  pairs 


w(r,s)" : 

(1.1)  ...  (0,  -sin<l>cos<}>,  r  cos2^  , 

(1.2)  ...  (cosatgij),  0,  0}  , 

(1.3)  ...  (cotga,  0,  0}  , 


\  (3.44) 

(2.2)  ...  {0,  (v2cos2a/r2)sin  26^,  (v2cos2a/r)cos26j)  , 

(2.3)  ...  (0,  (v2/r2)cosacotgacos  26j,  -(v2/r)cosacotgasin61cos61>, 

(3.3)  ...  (0,  -(v2/r2)cotg2asin  2^,  (v2/r)cotg2asin26j}  . 


The  results  (3.44)  are  to  be  contracted  with  aPPearin9 

in  (3.19).  The  thus  obtained  expressions  are  now  listed  according  to  the  (sym¬ 
metric)  pairs  "(i,j)":  'l 


(1.1)  ...  {0,  -tg^/r2,  1/r)  , 

(1.2)  ...  {sin((j>+61+x)/(r2cos2<()),  0,  0}  , 

(1.3)  ...  {-cos(^+61+T)/(r2cos2<{)),  0,  0}, 

(2.2)  ...  (0,  sin  2((5j+x)/r2,  cosz(61+x)/r)  , 

(2.3)  ...  (0,  -cos  2{61+x)/r2,  ^sin  2(61+x)/r}  , 

(3.3)  ...  (0,  -sin  2(s1+x)/r2,  sin2(61+x)/r}  . 


(3.45) 


In  the  second  and  third  expressions  of  (3.45),  <j>  has  been  replaced  by  ^+6^. 
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Finally,  the  second  partial  derivatives  in  the  first  term  on 

the  right-hand  side  of  (3.38)  are  contracted  with  fr_  £s_  appearing  in 

(i)  (j)  . 

(3.43)  and  the  elements  of  (3.45)  are  contracted  with  Fk  =  3F/3x  .  Thus 
the  desired  second-order  derivatives  of  F  expressed  in  spherical  coor¬ 
dinates  are: 


32F/3S(j)3S(j-)  =  [l/(r2cos2<j>)]32F/3A2-  ( tg<p/ r 2 )  3F/ 34>  +  (l/r)3F/3r, 

32F/3S(Y)3S(2)  =  [cos(61+T)/(r2cos'^)]32F/3A3$ 

-  [sin(61+r)/(r  cos$)]  3zF/3X3r 

+  [sin(^+6j+T)/(r2cos2^)]3F/3X  , 

32F/3s^Y)3s(3"j  =  [s-*n( 6i+t)4( r2COS^)3  32F/3X3^ 

♦  [cos(6j+T)/(r  coslj>)] 3zF/3X3r 
+  [cos(^+6^+x)/(r2cos2$)]3F/3X  , 

32F/3S(2)3s^2)  =  [cos2(61+T)/r2]32F/3$2  -  [sin  2(S1+x)/r]3zF/3?3r 
+  sin2(61+x)32F/3rz  +  [sin  Zfa^+xJ/r2] 3F/3^> 

+  [cos2(6j+x)/r]3F/3r  , 

32F/3S(2)3S(3)  =  [%sin  2(S1+x)/rz]3zF/3$2 +  [cos  2(<S1+x)/r]3zF/3$3r 

-  issin  2(61+x)3zF/3r2-  [cos  2(61+x)/r2]3F/3? 

+  *s[sin  2(61+i)/r] 3F/3r  , 


(3.46) 
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a2F/3s^3s^2)  =  [sin2(61+T)/r2]d2F/d<pz  +  [sin  2(6j+T)/r]32F/3$3r 
+  cos2(61+t)32F/3r2  -  [sin  2(61+T)/r2]3F/3$ 

+  [sin2(61+T)/r]3F/3r  . 

As  in  the  preceding  section,  if  t  is  set  to  62  ,  0  ,  or  -6j  the  desired 
derivatives  are  obtained  for  the  triads  i"  ,  i'  ,  or  i  ,  respectively. 
These  results  could  again  be  refined  upon  transformation  from  the  triad 
i"  to  the  triad  i  in  the  manner  of  equation  (3.34)  in  which  the  symbol 
F'  would  now  be  replaced  by  F. 
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tr  =  lr 

(T)  (1) 


{l/(r  cos<j»),  0,0}, 


£j_)  =  cos^j+t)^-  sin(61+T)£ 


(2) 


'(3) 


{0  ,  cos(61+t}/r, 
sin(61+T}}  , 


)  (3.49) 


lrm  =  sin(6.+T)e^  +cos(6,+t)£r  =  {0  ,  sin(6,+T)/r, 

(I)  1  (2)  1  (3)  1 

cos(6j+t)}  , 


which  agrees  with  (3.43). 

Upon  inserting  the  expressions  (3.47)  into  (3.20),  the  following 
nonzero  Christoffel  symbols  are  readily  obtained: 


=  » 


ri3  =  1/r  i 


’2  -  ei  n/f, 


11 


r2 

*23 


si n4>  cos<f>  , 


=  1/r  ; 


=  -r  cos*$  , 


r 00  —  ~  r  . 


22 


J 


(3.50) 
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The  application  of  (3.22)  yields  immediately  the  second  covariant  deriva¬ 
tives: 


Fji  =  32F/3A2-sin*cos*  3F/3?"+r  cos2^  3F/3r  , 
F12  =  32F/3A3?+tg$  9F/3A  , 

F13  =  a2,:/3X3r- (l/r>3F/3X  , 

F22  ■  32F/3<p  +  r  3F/3r  , 

F23  =  32F/3$3r-  (l/r)3F/3$  , 

F33  =  32F/3r2  . 


(3.51) 


With  the  aid  of  (3.49)  and  (3.51)  the  desired  second-order  derivatives 
32F/3s^-j3sqj  can  be  computed  according  to  (3.9).  Upon  carrying  out  the 
indicated  operations  the  results  (3.46)  are  obtained.  The  link  between 
the  two  coordinate  systems  established  in  the  preceding  section  has  thus 
confirmed  much  of  the  development  in  either  coordinate  system. 

By  setting  t  to  i2  ,  0  ,  or  -6j  the  desired  derivatives  are  ob¬ 
tained  for  the  triads  i"  ,  i'  ,  or  i  ,  respectively,  as  stated  earlier.  In 
order  to  evaluate  the  present  methodology,  the  case  x=-61  is  of  interest. 
Accordingly,  from  (3.46)  we  deduce 

32F/3S(1)3S(1)  =  [l/(r2COS^)]32F/3A2 

-  (tg^/r2)3F/3^+ (l/r)3F/3r  , 

32F/3S(i)3S(2)  =  [l/(r2COS$)]32F/3A3$ 

+  [sin<j>/(r2cos2$)]  3F/3A  , 


(3.52) 


32F/3s^1 j9s^3j  =  [l/(r  COS^)]  32F/3X3r  -  [l/(r2cos^)]  3F/3X  , 
32F/3s^2)3s^2)  =  (l/r2)32F/3$2+  (l/r)3F/3r  ,  / 

32F/3Sj2^3s^3j  =  (l/r)32F/3$3r  -  (l/r2)3F/3<jjT, 

32F/3s^3^3s^3j  =  32F/3r2  . 

Clearly,  this  outcome  could  be  obtained  also  directly  upon  using  (3.48)  and 
(3.51).  In  this  case,  only  the  steps  corresponding  to  (3.47),  (3.48), 
(3.50)  and  (3.51)  would  have  been  involved  while  (3.49)  and  the  more  gen¬ 
eral  results  (3.46)  would  have  been  by-passed. 

The  process  leading  to  (3.52)  would  have  required  considerably 
more  space  if  the  tensor  approach  to  this  problem  had  been  avoided,  whether 
for  pedagogical  or  other  reasons.  This  can  be  appreciated  upon  consulting 
the  development  in  [Tscherning,  1976],  where  our  formulas  (3.52)  were  ob¬ 
tained,  in  the  present  order,  in  equations  numbered  as  (53),  (55),  (57), 
(52),  (56)  and  (54);  our  notations  ds^  ,  ds^  *  ds(3)  corresPond  re¬ 
spectively  to  dy  ,  dx  ,  dz  in  this  reference.  The  methodology  used  there¬ 
in  was  not  based  on  tensor  analysis  but  made  use  of  the  ordinary  differen¬ 
tiation  and  interrelationships  between  the  local  Cartesian  coordinate  sys¬ 
tem,  the  earth-fixed  Cartesian  coordinate  system  and  the  polar  coordinate 
system.  Such  a  paper  is  certainly  useful,  especially  from  the  pedagogical 
point  of  view,  since  it  addresses  itself  to  large  audiences  and  exposes 
the  treated  subject  in  great  depth.  On  the  other  hand,  in  [Tscherning, 
1976]  the  derivations  leading  to  our  formula  (3.52)  required  eight  pages 
(pages  73  through  80),  even  though  only  the  simplest  case  corresponding 
to  the  triad  1  was  treated.  We  might  thus  conclude  that  although  the 
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methodology  of  Hotine  [1969]  requires  prior  knowledge  of  certain  elements 
of  tensor  analysis,  this  disadvantage  (if  it  is  indeed  a  disadvantage)  is 
more  than  offset  by  the  ease  and  rapidity  with  which  many  general,  as  well 
as  specialized,  formulas  in  geodesy  may  be  derived. 

The  well-known  formula  giving  the  Laplacian  of  F  in  spherical 
coordinates  is  readily  obtained  as  the  trace  of  the  matrix  composed  of 
the  elements  presented  in  (3.46)  for  any  choice  of  r;  the  easiest  choice 
is,  of  course,  that  of  t  = -<Sj  corresponding  to  (3.52).  In  every  case 
we  obtain 

AF  =  [l/(r2cos2$)]a2F/aX2+ (l/r2)32F/a^  +  32F/3r2 

-  (tg$’/r2)3F/3$'+  (2/r)3F/3r  .  (3.53) 

As  has  already  been  pointed  out  in  the  closing  statement  of  the  preceding 
section,  the  results  (3.46),  specialized  for  the  triad  i"  ,  can  be  further 
refined  upon  transformation  to  the  triad  i  ;  it  is  immediately  verified 
(see  3.34)  that  AF  remains  unchanged. 

A  practical  question  left  unanswered  so  far  concerns  the  deter¬ 
mination  of  the  angle  needed  in  specializing  the  results  (3.46)  for 
the  triad  iu  .  In  fact,  the  angle  needed  in  this  task  is  the  total  angle 

S-Aj+fig*  (3. 54) 

where  6j  is  known  from  (3.16).  In  order  to  determine  6  we  proceed  as  in 
Section  2  where  alone  was  sought.  The  vector  Up  will  now  be  projected 
onto  the  triad  (3.48),  the  projections  being  respectively 


1 


3U/3s^jj  =  =  0  , 

3U/3S(2)  =  Ur£[2)  =  (l/r)3U/34> 
3ll/3s^2 j  =  =  • 


(3.55) 


Similar  to  (3.29)  we  have 


tg6  =  (-3U/3s,9.)/(-3U/3sm)  , 


(3.56) 


which  is  in  general  positive  whether  in  a  rotating  or  a  nonrotating  field 
(at  the  poles  and  on  the  equator  it  is  zero).  The  radial  (positive  in¬ 
ward)  and  the  meridian  (positive  north)  components  of  Y  at  P  are  now  de¬ 


noted  as  yp  ,  Ym  ;  we  have 


*n  "  -3U/3s(3)  . 


Ym  -  3U/3s(2)  , 

Y  =  (Y*  +  Y*)55  , 
n  m 


and  from  (3.56), 


tgs  -  -Ym/Yn 


(3.57a) 

(3.57b) 


(3.57c) 


(3.58) 


In  expressing  (3.55)  explicitly,  use  can  be  made  of  the  spherical 
harmonic  expansion  of  U  which  reads  (see  e.g.  [Blaha,  1977],  page  79) 


U  =  (GM/r)[l  +  (ao/r)2C*0P2(sin?)  +  (ao/r)**C|0P4(sin?)  +  ...] 


+  *s&2r2cosz<p  , 


(3.59) 


where  C|g  ,  C|g  ....  ,  are  the  spherical  harmonic  coefficients  in  the 
standard  field,  linked  to  iGA2g  in  (3.27b)  as  follows: 

» 

C|2n-2),0  =  (-Dn[-(2n  +  l)  +  (2n-2l  Q  e2n_2  ]  /  [(2n  -  l)(2n  +  1)]  , 

Q  =  iGA20(ae)/GM  . 

In  analogy  to  Section  2,  if  one  wishes  to  consider  a  nonrotating  field  the 
terms  containing  £2  explicitly  should  be  dropped.  From  (3.59)  we  express 

3U/3*  =  (GM/r)  [(aQ/r) 2C*gdP^sih$)/#+  (a0/r)4C|gdP4(sin?)/#+  . . .] 

-  S2r2siri^cos^  ,  (3.60) 

3U/3r  =  -(GM/r2)  [l  +  3(aQ/r) 2C|QP2(sin?)  +  5(ao/r)4C*gP4(sin^)  +  . . .] 

+  5i2r  cos2^  ;  (3.61) 

the  terms  containing  the  coefficients  beyond  C|g  and  certainly  those  beyond 
C£q  are  usually  neglected.  For  reference  purposes  we  list: 

P2(sin$)  =  *s(3  sin2$-  1)  , 

P4(sin$)  =  (1/8) (35  sin4^-30  sin2^+3)  ; 

dP2(sin$)/d<£=  3  sin$cos$  , 
dP4(sin<j>)/d$  =  ( 5/ 2 ) ( 7  sin2$- 3)sin^cos^  ; 

d2P2(siri$)/d"ip  =  3(-2  sin2$+ 1)  , 
d2P4(sin$)/dF  =  Js(-140  sin4$+135  sin2*- 15)  . 

The  inclusion  of  C?n  would  follow  along  the  same  lines. 


J.S  l’rd(.J.1c.j  I  AjjjjI  ii.ij|)l  I  I  Ly  of  l<i“.ull.s 

In  this  section  the  function  F  is  treated  in  spherical  coordin¬ 
ates  and  is  first  specialized  to  U  (standard  potential)  and  then  to  T 
(disturbing  potential).  Also,  the  more  usual  geodetic  conventions  cor¬ 
responding  to  (3.3)  are  adopted.  First  we  consider 

F  =  U  , 

r  =  S2  (the  triad  i  is  replaced  by  the  triad  i")  ; 

taking  into  account  the  property  9U/9A=0  and  the  convention  (3.54),  the 
results  (3.46)  read 

U  „  „  =  (cos26/r2)92U/9?2  -  (sin  26/r)92U/9$9r 

A  A 

+  sinz6  32U/9r2  +  (sin  26/r2)dU/9$  +  (cosz6/r)dU/dr  , 

Ux"y"  =  0  ’ 

Ux„z„  =  (h  sin  26/r2)92U/9$2 +  (cos  26/r)92U/9^9r 

-  h  sin  26  92U/9r2  -  (cos  25/rz)9U/9$ + %(sin  26/r)9U/9r  , 

Uyy,  =  -(tg4>/rz)9U/94>+  (l/r)9U/9r  ,  >  (3.63) 

*V'z"  =  ^  ’ 

Uz„z.i  =  (sin26/r2)  92U/9^2  +  (sin  26/r)9zU/9$9r 

+  cosz6  92U/9r2-(sin  26/rz)9U/94>+ (sin26/r)9U/9r  . 
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Equation  (3.53)  applies  now  in  the  form 


AU  =  Ux„x„  +  Uy„y(l  +  Uz„2„  =  (l/r2)32U/3?2  +  32U/3r2 

-  (tg?»/r2)9U/3?+ (2/r)3U/3r  .  (3.64) 

The  first  partial  derivatives  needed  in  (3.63),  (3.64)  have  been  presented 
in  (3.60)  -  (3.62).  The  second  partial  derivatives  follow  as 


32U/3^  =  (GM/r)  [(aQ/r) 2C|0d2P2(sin<p)/d«p2 

+  (a0/r)ltC|od2P4(sin^')/d^2  +  . . .]  -  a2r2(cos2<jT-  sin2<jO  , 


32U/3$3r  =-(GM/r2)[3(ao/r)2C*0dP(sin$}/d$‘  S  (3.65) 

+  5(ao/r)4C|0dP4(sin$)/d$'+ . . .]  -  2532r  sin$cos<jr  , 

32U/3r2 =  2(GM/r3) [1 +  6(aQ/r)2C|0P2(sin^) 

+  15(ao/r)4C|0P4(sin?))  +  . . .]  +  S2cos2?  , 


in  which  the  formulas  (3.62)  apply  again.  If  these  explicit  partial  deriva¬ 
tives  of  U  are  used  in  (3.64),  the  coefficients  of  C£Q  ,  ,  ...  ,  are 

seen  to  be  zero  and  the  Laplacian  is  confirmed  to  be 

AU  =  2a2  .  (3.66) 

If  a  nonrotating  field  is  considered  the  Laplacian  becomes  zero  and  U  be¬ 
comes  a  harmonic  function  (as  is  T  or  the  gravitational  part  of  W). 


If  we  wish  to  obtain  the  second-order  derivatives  of  U  in  the 


triad  i  rather  than  in  the  triad  i"  ,  the  expressions  in  (3.34)  are 
specialized  (U  is  written  for  F'  ,  etc.)  and  combined  with  (3.63)  to 
give 

(3.67) 

Uxx  ~  Ux“x"  "  Ux"z"  =  [(cos26-C  sin  26)/r2]32U/34>2 
-[(sin  26  +  2£cos  26)/r]  32U/3?3r 

+  (sin2S  +  £sin  26)32U/3r2 +  [(sin  26  +  2Ccos  26)/rz]3U/3^ 

+  [(cos26-£sin  26)/r]3U/3r  , 

Uxy  “  "nUx"z"  =  ‘nWs1n  26/r2)32U/3?2  +  (cos  26/r)32U/3^3r 

-  ^  sin  26  32U/3r2  -  (cos  26/r2)3U/3^  + ^(sin  26/r)3U/3r]  , 

Uxz  =  Ux"z"  ‘  ?(Uz"z"  "  Ux"x")  =  Sln  26  +  ^cos  26)/r2]32U/3^2 
+  [(cos  26  -  2£sin  26)/r  ]  32U/3?3r  -  sin  26  +  £cos  26)3zU/3r2 

-  [(cos  26-2?sin  26)/r2]  3U/3^  +  [(^sin  26  +  £cos  26)/r]3U/3r  , 

Uyy  2  Uy"y"  =  -(tg4>/r2)3U/3l£+  (l/r)3U/3r  , 

U~2  -  -n(Uzllz„  -  U  „  =  -n((sin26/r2)32U/3^2  +  (sin  26/r)32U/3$3r 

+  C0S26  32U/3r2+  [(tgif-  sin  26)/rz]3U/3?  -  (cos26/r)3U/3r>  , 

U22  ~  Uz"z"  +  Ux"z"  =  [(Sin2<s  + ^sin  26)/r2]32U/3^2 

+  [(sin  26 + 2£cos  26)/r]  32U/3^3r +  (cos26  -  ^sin  26)3zU/3r2 
-[(sin  26  +  2Ccos  26)/r2]3U/3$  +  [(sin26  +  £sin  26)/r]3U/3r  . 
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For  practical  applications,  the  relation 


W  =  U  +  T 


leads  to  =  U~|  +  ,  etc.,  where  W~~  ,  etc.,  may  be  directly  observed. 

The  observations  could  be  considered  at  the  ground  level  or  at  higher 
levels  (aircraft  altitude,  satellite  altitude)  and  distinction  could  be 
made  between  the  rotating  and  nonrotating  fields.  From  the  observations 
W-j  »  etc.,  and  from  their  counterparts  ,  etc.,  the  values  »  etc., 
can  be  derived  and  modeled  by  spherical  harmonics,  point  masses  and  other 
means.  For  example,  lower  degree  and  order  spherical  harmonic  coefficients 
may  be  considered  fixed  and  higher  degree  and  order  coefficients  may  be 
solved  for  in  a  least  squares  adjustment;  or  instead  of  solving  for  these 
coefficients,  local  parameters  such  as  point  mass  magnitudes  can  be  adjust¬ 
ed.  In  the  most  rigorous  case,  the  values  to  be  modeled  are 


XX 

XX 

XX 

T~~  = 

W~~  - 

u~~ 

xy 

xy 

xy 

= 

- 

xz 

xz 

xz 

T~~  = 

- 

U~~ 

yy 

yy 

yy 

T~~  * 

w~~  - 

u~~ 

yz 

yz 

yz 

Tz2  * 

Wzz  ■ 

Uzz 

(3.68) 
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Mni  u  I  ,  <in*  n*l#il.  Iv<*ly  *.m#i  1 1  guiinl  1 1  Ic. ,  lh«*  triad  I 

involved  In  their  formulation  tan  be  replaced  by  the  triad  I"  ,  i‘  , 
or  even  i  ;  the  case  i'  can  be  viewed  as  a  spheroidal  approximation  and 
the  case  i  ,  as  a  spherical  approximation.  On  the  other  hand,  the 
spherical  approximation  with  regard  to  ,  etc.,  would  in  general  re¬ 
sult  in  prohibitive  errors.  If  a  simplification  is  to  take  place  in 
this  respect,  we  assume  that  it  will  result  in  Uz„z„  ,  etc.,  replacing 
>  etc.  Below  are  listed  the  practical  counterparts  of  (3.68).  With 
regard  to  T,  the  triad  i  identifying  the  spherical  approximation  (in  paren¬ 
theses)  and  the  triad  i"  are  considered: 


(Txx  ^  Tx"x" 
( Tx"y" 
(Txz^  Tx"z" 


(Tyy  ^  Ty"y" 
^Tyz~)  Ty"z" 
<Tzz=>  Tz„z„ 


W~-w  —  U  it  ii  i 

XX  XX 


~  w 

xy  ’ 


=  "  UxMz"  ’ 


Wyy  "  Uy"y"  ’ 


Wyz  ’ 


-  Wzz  -  Uz"z"  ’ 


•\ 


(3.69) 


the  values  Uzl,z„  ,  etc.,  have  been  presented  in  (3.63)  and  the  values 

Tzz  *  Tz"z"  ’  etc”  wil1  be  develoPed  shortly.  It  can  be  shown  that 

errors  caused  by  simplifications  in  the  horizontal  gradients  of  standard 

gravity  (U~~  ,  U-~)  and  in  the  vertical  gradient  of  standard  gravity  (U~~) 
yz  z z 

reach  approximately  the  following  values: 
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0.011  E  per  1"  in  £  , 

0.011  E  per  1"  in  n  » 

0.000024  sin  2$  E  per  1"  in  £  . 

The  error  in  U~~  is  completely  negligible  under  any  circumstances  since 
even  if  £  were  60"  it  would  only  reach  0.001  E  (the  measurements  are 
often  considered  good  to  ±0.05  E,  where  E e  Eotvos  =  0.1  mgal/km). 

In  order  to  develop  the  second-order  derivatives  of  T  in  (3.69), 
we  specialize  (3.46): 

F  =  T  , 

T  E  62  • 

The  expressions  for  Tz„z„  ,  etc.,  can  be  obtained  by  simply  transcribing 
(3.46)  with  T  replacing  F  ,  6  replacing  (6j+t)  and  dx,  etc.,  replacing 
ds^)  »  etc.;  there  is  no  need  to  list  the  results  again  with  these  minor 
notational  changes.  When  T  is  developed  in  spherical  harmonics,  such  re¬ 
sults  could  be  made  computationally  more  economical  by  eliminating  a2T/9^2 
which  involves  the  evaluation  of  d2Pnm(sin^')/d^2  ,  where  Pn[n(sin<f)  are  the 
associated  Legendre  functions  in  argument  sinf.  By  contrast,  when  U  was 
considered,  the  corresponding  two  or  at  most  three  second  partial  deriva¬ 
tives  with  respect  to  ^  were  exceedingly  simple  to  evaluate  in  terms  of 
even  powers  of  sin<£  (see  equation  3.62).  In  the  present  situation,  a2T/a^2 
can  be  eliminated  upon  using 

AT  =  0,  (3.70) 


Jx2  ”• 


yz 


22 
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wher*'  the  Laplacian  of  T  is  given  in  (3.53)  with  T  replacing  F  ,  hence 


( \/'^'  -  -  Ll/fr^cos^JJa^T/aA2  -  j'T/ar2 

+  (tg^/r2)3T/3$  -  (2/r)3T/3r  .  (3.71) 


With  the  substitution  (3.71)  the  transcription  of  (3.46)  yields 
upon  incorporating  the  above-mentioned  notational  changes: 

Tx.x-  =  -  [cos26/(r2cos2$)]32T/3A2  -  (sin  26/r)32T/3lj>3r  " 

-  cos  26  32T/3r2  +  [(tg$  cos26  +  sin  26)/r2]3T/3l£ 

-  (cos26/r)3T/3r  , 


T  „  „  =  [cos6/(r2cos<p)J  32T/3X3<f>  -  [sin6/(r  cos<)>)]  32T/3X3r 
*  y 

+  [sin(‘$+6)/(r2cos2^)]  3T/3X  , 

Tx,.z„  =  -%[sin  26/(r2cos2?)] 32T/3X2  +  (cos  26/r)32T/3^3r 

-  sin  26  32T/3r2  -  [(cos  26-%  sin  26  tg^T)/r2]3T/3^ 

-  %(sin  26/r)3T/3r  , 


(3.72) 


Vy"  =  f1/^20052^)]921/^2-  (tg^/r2)3T/3?+(l/r)3T/3r  , 

TyMz„  =  [sin6/(r2cos^)]32T/3X3^'  +  [cos6/(r  cos^)]32T/3X3r 

-  [cos($+6)/(r2cos2^)]3T/3X  , 

Tz"z"  =  -[sin26/(r2cos2^)]32T/3X2 +  (sin  26/r)32T/3$3r 
+  cos  26  32T/3r2  -  [{sin  26  -  sin26tg^)/r2]3T/3? 

-  (sin26/r)3T/3r  . 


Some  of  the  algebra  is  verified  immediately  by  confirming  (3.70). 
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The  specialization  of  (3.46)  with  t  h  -fij  already  resulted  in  (3.52). 
If  also  (3.71)  is  utilized  (3.52)  can  be  transcribed  as 


Txx  =  .  [i/(r2cos2$)]32T/ax2  -  a2T/ar2  +  (tg<f>/r2)aT/a* 

-  (l/r)3T/ar  , 

T  =  [l/(r2cos?)3a2T/axa?+ [tg?/(r2cos^)]aT/aX  , 
xy 

Txz  =  (l/r)a2T/ai$ar-  (l/r2)aT/a*  ,  > 

T  =  [l/(r2cos2^)]a2T/ax2  -  (tg$/r2)aT/a$ +  (l/r)3T/ar  , 

yy 

Tyz  =  cos?)]a2T/axar-  [l/(r2cos?)]3T/aX  , 

Tzz  =  a2T/ar2  .  , 


(3.73) 


However,  in  this  case  only  Txx  originally  contained  a2T/a<|>2  ;  this  original 
form  is 

Txx  *  d/r2)a2T/a?2+ (l/r)3T/ar  .  (3.73') 


Before  discussing  the  usefulness  of  ( 3 . 73 1 )  and  some  of  the  circum¬ 
stances  leading  to  the  choice  of  the  triad  i  over  the  triad  1“  or  vice  versa, 
we  present  yet  another  set  of  formulas  which  may  be  said  to  characterize  an 
intermediate  precision.  This  set  is  essentially  (3.72)  where  the  second 
and  higher  powers  of  6  are  neglected  in  the  series  expansion  of  trigonometric 
functions.  We  thus  obtain  (with  3.71  already  incorporated)  from  either 
(3.72)  or  (3.73): 
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Vx"  *  Txx~26Txz  =  -[l/(r2C0S2$)]92T/3A2- 2(6/r)32T/3$3r 
-  32T/3r2+ [(tg$  +  26)/r2]3T/3$- (l/r)3T/3r  , 

Tx"y"  ~  Txy“6Tyz  =  [V(r2cos^)]32T/3A3$'-  [6/ ( r  cos?)]  32T/3X3r 
+  [(sin4>  +  5cos$)/(r2cos2$)]  3T/3A  , 


T 


x"z' 


*  Txz  +  6(Txx-T22)  =  -[6/(r2cos2^)]32T/3A2+ (l/r)32T/3^3r 
-  26  32T/3r2  -  [(1  -  6tg$)/r2] 3T/3$  -  (6/r)3T/3r  , 


>(3.74) 


Vy*  =  Tyy  =  [l/(^2cos2^)]  32T/^X2  -  (tg<j>/r2)3T/3<j>+  (l/r)3T/3r  , 

V'z"  ~  Tyz  +  5Txy=  [,$/(r2cos^)]  32T/3A3<F+  [l/(r  cos<jT)]  32T/3A3r 
-  [(cosiji'- 6sin$)/(r2cos2<F)]3T/3A  , 


Tz"z"  “  Tzz  +  25Txz  =  2(6/r)32T/3<j>3r  +  32T/3r2  -  2(6/r2)3T/3(jT  .  J 


The  angle  6  called  for  in  (3.72)  or  (3.74)  is  usually  known  beforehand  since 
it  is  needed  for  computing  the  second-order  derivatives  of  1)  in  (3.63).  If 
this  angle  were  not  readily  available  (e.g.  if  the  second-order  derivatives 
of  T  were  supplied  directly  and  would  not  have  to  be  determined  from  the 
second-order  derivatives  of  W),  it  could  be  approximated  by  6j  =  <J>-4>  and 
we  would  then  be  in  the  presence  of  the  spheroidal  approximation  mentioned 
earlier.  In  this  case  one  would  quite  logically  use  (3.74)  over  (3.72), 
since  the  errors  caused  by  the  spheroidal  approximation  would  be  in  general 
much  larger  than  the  errors  in  (3.74)  caused  by  the  simplification  of  the 
formulas  in  (3.72). 


It  has  been  suggested  that  the  quantities  to  be  modeled  (second- 
order  derivatives  of  T)  could  be  described  by  more  than  one  set  of  para¬ 
meters,  for  example  by  a  given  set  of  spherical  harmonic  coefficients  and 
a  set  of  local  parameters  such  as  point  mass  magnitudes.  The  first  set 
could  be  considered  known  and  thus  only  the  second  set  (e.g.  point  mass 
magnitudes)  could  be  subject  to  adjustment.  This  second  set  would  be  re¬ 
quired  to  accommodate,  in  a  least  squares  adjustment,  the  "residual  second- 
order  derivatives"  of  T  obtained  by  subtracting  the  "computed  second-order 
derivatives"  of  T  from  the  original  ones.  The  "computed"  quantities,  ob¬ 
tained  via  fixed  spherical  harmonic  coefficients,  could  be  evaluated  as  in 
(3.72);  the  partial  derivatives  32T/3A2  ,  etc.,  will  be  listed  below.  The 
"computed"  quantities  can  be  easily  obtained  without  introducing  any  errors, 
thus  (3.72)  is  indeed  suitable  for  this  purpose.  The  formulas  (3.74)  of 
intermediate  precision  would  be  almost  equally  suitable.  It  is  conceivable 
that  for  many  practical  applications  even  the  formulas  (3.73)  could  be  ac¬ 
ceptable.  On  the  other  hand,  when  dealing  with  the  "residual  second-order 
derivatives"  the  formulas  (3.73)  are  believed  to  be  sufficient  in  most  cir¬ 
cumstances  and  could  be  used  with  great  advantage.  For  example,  in  the  case 
of  a  point  mass  model  important  mathematical  simpl ications  will  occur  if 
these  "residual"  quantities  are  described  by  (3.73)  provided  they  can  be 
considered  distributed  on  or  near  a  sphere  and  the  point  masses  are  chosen 
to  lie  on  another  (concentric)  sphere  of  a  smaller  radius.  It  may  thus  be 
necessary  to  introduce  further  simplifications  in  the  formulas  (3.73)  in 
order  to  make  them  computationally  viable  for  an  actual  adjustment  process. 
One  notices  that  in  the  case  of  point  masses  or  other  local  parameter, 
there  is  no  reason  to  avoid  the  expression  3*17  3$2  and  it  may  be  more  ad¬ 
vantageous  to  use  (3.73')  in  lieu  of  its  counterpart  in  (3.73). 


For  reference  purposes,  the  partial  derivatives  of  T  with  re¬ 
spect  to  the  spherical  coordinates  are  now  listed  for  T  expanded  in  terms 
of  spherical  harmonics.  As  explained  earlier,  92T/9$2  is  not  among 
them.  The  familiar  expression  for  the  disturbing  potential  is  (see  e.g. 
[Blaha,  1977],  page  80) 


N 

T  =  (GM/r)  l  (a  /r)n  AS(n)  , 
n=2  u 

where 

n 

AS(n)  =  l  (AC  cos  mX  +  AS  sin  mx)P  m(siri$)  . 
m=0 


(3.75) 


(3.76) 


N  =  degree  and  order  of  the  truncation  , 


AC20  =  C20  '  C20  ’ 

AC40  =  C4Q  “  C4Q  »  etc- » 


the  other  AC  ,  AS„„  being  C  ,  S  m  ,  the  spherical  harmonic  potential 
nm  nm  3  nm  nm 


coefficients  themselves  (in  practice  this  could  apply  already  for  C6Q). 
We  then  have 


9T/9X 


N 

(GM/r)  l  (a/r)n  AS’(n)  , 
n=2  0 

N 


9T/9$  =  (GM/r)  l  (an/r)n  AS(n)  , 


n=2 

N 


9T/9r  =  -(GM/r2)  £  (n + 1 )(a  /r)n  AS(n)  , 

n=2 

N 

92T/9X2  =  (GM/r)  £  (a  /r)n  AS"(n)  , 

n=2  0 
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(3.77) 


32T/3X3*  * 
32T/3X3r  = 
32T/3?3r  = 
32T/3r2  = 


(GM/r)  1  (a Q/r)n  AS ' ( n )  , 
n=2 

N 

■(GM/r2)  l  (n  +  l)(a  /r)n  AS ‘ ( n )  , 
n=2  0 


N 


■(GM/r2)  l  (n  +  l)(a  /r)n  AS(n)  , 
n=2  0 

N 

(GM/r3)  l  (n+ l)(n  + 2)(a  /r)n  AS(n)  , 
n=2  0 


where  the  following  four  notations  have  been  introduced: 

n 


AS' (n) 
AS"(n) 
AS(n) 
AS' (n) 


=  m=0  m('4C<»Sin  mX  +  4S™c0s  "rt)pral(si"«  • 
n 

*  -7  m2(AC  xos  mX  +  AS  sin  mX)P  (sin4>)  , 

m^Q  nm  nm  nm 

n 

*  Jo  (ACnmcos  mX  +  AS^sin  mX)dPnm(sin^)/#  , 

n  _ 

=  T  m(-ACsin  mX  +  AS„mcos  mX)dP  (sin<f>)/d<j>  . 
*•„  nm  nm  nm  T 

m=u 


(3.78a) 
(3.78b) 
(3.78c) 
(3. 78d) 


3.6  Summary  and  Conclusion 


One  vertical  and  two  horizontal  gradients  of  gravity  are  among 
the  quantities  which  can  routinely  be  measured  and  used  in  specific  geo¬ 
detic  applications.  They  represent  three  out  of  six  distinct  second-order 
derivatives  of  the  potential  in  a  local  Cartesian  coordinate  system  as  re¬ 
vealed  by  various  sensors.  These  six  quantities  describe  the  curvature 
properties  of  the  gravity  field.  If  their  knowledge  is  to  contribute  to 
a  better  description  of  the  earth's  gravity  field,  they  should  be  properly 
related  to  the  model  parameters  of  the  field,  such  as  a  set  of  spherical 
harmonic  potential  coefficients  and/or  a  set  of  chosen  localized  parameters, 
etc.  However,  before  these  second-order  derivatives  can  be  related  to  a 
set  of  parameters,  they  have  to  be  expressed  in  terms  of  the  partial  deri¬ 
vatives  of  the  potential  with  respect  to  the  coordinates  of  a  chosen  coor¬ 
dinate  system.  In  this  chapter,  the  coordinate  systems  used  are  the  spheroi¬ 
dal  coordinate  system  {<$,u,a}  and  the  spherical  coordinate  system  {X,f,r}. 

The  model  parameters  themselves  are  brought  into  the  picture  only 
in  the  closing  paragraph  of  Section  3.5  where  they  are  limited  to  a  set  of 
spherical  harmonic  potential  coefficients.  For  the  most  part,  this  chapter 
has  been  concerned  with  the  second-order  derivatives  of  a  general  scalar 
function  of  position  (F)  with  respect  to  the  length  elements  along  a  family 
of  local  Cartesian  axes,  expressed  in  terms  of  the  partial  derivatives  of  F 
with  respect  to  the  coordinates  of  the  two  systems  just  mentioned.  The  func¬ 
tion  F  is  most  likely  to  represent  W  (actual  potential),  U  (standard  poten¬ 
tial),  T  (disturbing  potential),  or  a  coordinate  itself  (a,  r,  etc.). 
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The  desired  second-order  derivatives  of  F  are  expressed  in  sphe¬ 
roidal  coordinates  in  Section  3.2  and  in  spherical  coordinates  in  Section 
3.4.  A  link  between  the  two  formulations  is  established  in  Section  3.3 
where  all  the  results  in  spherical  coordinates  are  obtained  also  indirectly 
through  a  transformation  from  the  spheroidal  system.  In  addition  to  pos¬ 
sible  applications  of  similar  links,  a  by-product  of  such  a  demonstration 
is  the  verification  of  the  results  obtained  in  either  system  by  direct 
methods. 

A  benefit  which  surfaced  during  the  development  in  spherical  coor¬ 
dinates  is  that  the  demonstrations  can  be  significantly  shortened  and 
brought  sharply  into  focus  if  the  approach  pioneered  by  Marussi  [1951]  and 
elaborated  by  Hotine  [1969]  is  followed  and  advantage  is  taken  of  tensor 
analysis.  It  thus  comes  as  no  surprise  that  this  study,  in  which  such  a 
methodology  has  been  followed,  is  intended  as  a  tribute  to  Professor 
Marussi  and  to  the  late  Martin  Hotine  for  their  contributions  to  theoretical 
geodesy. 

The  general  formulas  giving  the  second-order  derivatives  of  F  are 
made  useful  for  practical  applications  in  Section  3.5,  upon  decomposing  W 
into  U  and  T.  At  least  some  of  the  second-order  derivatives  of  W  are  con¬ 
sidered  as  observed  quantities  and  F  is  first  specialized  to  represent  U. 

In  this  way  the  second-order  derivatives  of  T =  W-U  are  obtained  and  may  in 
turn  be  considered  as  observations  to  be  modeled  by  a  chosen  set  of  para¬ 
meters.  For  this  purpose  the  new  "observations"  are  expressed  in  terms  of 
the  partial  derivatives  of  T  with  respect  to  the  coordinates,  accomplished 
by  specializing  F  to  represent  T.  The  above  practical  procedure  is  developed 
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in  detail  in  spherical  coordinates.  If  the  computation  of  32T/3<p  is  more 
complicated  than  that  of  the  other  partial  derivatives  --  as  would  be  the 
case  with  the  spherical  harmonic  expansion  —  this  computation  can  be  avoid¬ 
ed  by  taking  advantage  of  the  property  AT =  0.  A  typical  set  of  second-order 
derivatives  of  T  with  the  property  AT  =  0  incorporated  (and  32T/3q>2  absent) 

can  be  found  in  (3.72).  If  there  is  no  need  to  avoid  the  formulation  of 

< 

32T/ 3<ji2,  the  desired  second-order  derivatives  can  be  readily  transcribed 
from  equations  (346). 

The  possibility  of  considering  a  whole  family  of  local  Cartesian 
coordinate  systems  has  lead  to  the  second-order  derivatives  of  T  in  the 
following  two  basic  forms.  The  first  form  (associated  with  the  Cartesian 
triad  called  i")  is  essentially  rigorous  and  is  expressed  in  the  above- 
mentioned  equations  (3.72)  or,  eventually,  (3.46).  The  second  form  (asso¬ 
ciated  with  the  Cartesian  triad  called  i)  corresponds  to  the  spherical  ap¬ 
proximation  and  is  expressed  in  equations  (3.73)  where  AT =  0  has  been  in¬ 
corporated;  if  there  is  no  need  to  avoid  the  formulation  of  32T/34>2,  the 
first  relation  in  (3.73)  can  be  replaced  by  (3.73').  Another  form  of  inter¬ 
mediate  precision  (with  AT =  0  incorporated)  may  be  found  in  (3.74). 

The  formulas  in  this  chapter  have  been  developed  for  all  six 
second-order  derivatives  of  W,  U,  T,  or  a  more  general  function,  although 
only  one  vertical  and/or  two  horizontal  gradients  of  gravity  may  be  routine¬ 
ly  measured.  One  notices  that  if  the  observed  quantities  were  given  di¬ 
rectly  as  the  second-order  derivatives  of  T  (instead  of  W),  the  approach 
and  all  the  formulas  would  remain  unchanged  except  that  the  steps  associated 
with  U  would  be  skipped. 
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The  versatility  in  possible  applications  of  the  presented  formu¬ 
las  becomes  further  apparent  upon  the  realization  that  they  have  been  deriv¬ 
ed  for  a  general  point  in  space  and  can  thus  be  used  for  measurements  made 
at  the  ground  level  as  well  as  at  higher  altitudes,  such  as  aircraft  or 
satellite  altitudes. 

Although  the  formulas  giving  the  desired  second-order  derivatives 
are  to  be  applied  most  often  in  a  rotating  field,  they  can  be  easily  spe¬ 
cialized  for  a 
applications, 
of  the  earth's 
citly,  and  the 
orientation  of 


nonrotating  field  as  could  be  required  for  certain  satellite 
This  specialization  consists  in  replacing  ai(angular  velocity 
rotation)  by  zero  wherever  it  appears  in  the  formulas  expli- 
same  modification  also  takes  place  in  connection  with  the 
the  triad  i". 
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'I.  AIjJIJ;>  IMLN  f  MOUl.Lfi  I  OR  I  IVL  INVLO  I  IOAI  I.L) 
OBSERVATIONAL  MODES 

4.1  Mathematical  Background 

In  a  simultaneous  adjustment  of  spherical  harmonic  (S.H.)  poten¬ 
tial  coefficients  and  point  mass  (P.M.)  magnitudes  as  parameters,  an  obser¬ 
vation  equation  would  read  as  follows: 


F(xJ,Xa)  -  Lb 

(4.1a) 

F(Xa ,0)  +  A2X2  -  Lb 

(4.1b) 

F<X°,0)  +  A1X1  +  A2X2  -  Lb  , 

(4.1c) 

where  V  is  the  residual,  La  and  Lb  are  the  adjusted  and  measured  values  of 
the  observable,  the  subscripts  "1"  and  "2"  refer  to  the  S.H.  and  P.M.  para¬ 
meters,  and  the  superscripts  "a",  "o"  attributed  to  "X"  indicate  the  adjusted 
and  initial  values  of  the  parameters.  Furthermore,  "F"  denotes  the  general 
model  describing  the  observable,  "A."  denotes  the  row  matrix  of  partial  deri¬ 
vatives  of  F  with  respect  to  the  appropriate  parameters  (sets  "1"  and  "2" 
above),  and 

Y  _  yO 

A1  “  A1  "  A1  > 

X2  =  x|  ,  x£  =  0  ; 

the  last  relation  indicates  that  the  model  is  linear  in  the  P.M.  parameters. 
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Equation  (4.1b)  is  written  as 


V  =  A2X2  +  [F(xj,0)-Lb]  ,  (4.2) 

where 

F(X*,0)  =  F(X° ,0)  +  AjXj  ; 

the  expression  in  brackets  of  (4.2)  becomes  the  constant  term  in  the  obser¬ 
vation  equation  in  the  cases  where  the  set  X®  is  given  (fixed)  and  only  X2 
is  subject  to  adjustment.  In  either  case,  the  expressions  to  be  developed 
include  F(X°,0)  and 

AFj  =  AjXj  , 

AF 2  A 2^2  * 

In  practice  it  is  useful  to  proceed  in  two  consecutive  adjustments 
of  the  same  observations.  First,  one  minimizes,  in  the  least-squares  (L.S.) 
sense,  the  expressions  in  the  brackets  of  (4.2)  as  if  the  observation  equa¬ 
tions  were  of  the  type 

Vj  =  AjXj  +  [F(X°,0)  -  Lb]  ,  (4.3a) 

and  then  one  forms 

V  =  A2X2  +  Vx  ,  (4.3b) 

as  if  Vj  represented  the  minus  observed  values  in  a  linear  model  where  X2  are 
the  parameters.  The  equations  (4.3a)  and  (4.3b)  considered  together  are 
identical  in  form  to  (4.1c).  However,  the  adjustment  itself  is  performed  in 


two  steps  rather  than  simultaneously  as  suggested  by  (4.1c).  Splitting 
the  adjustment  in  this  way  is  possible  in  some  cases  without  introducing 
any  approximations;  for  example,  with  the  S.H.  coefficients  as  parameters 
and  with  an  almost  perfect  distribution  of  observations  (such  as  geoid  un¬ 
dulations)  over  the  entire  globe,  two  or  more  subsets  of  these  parameters 
can  be  solved  for  consecutively  with  no  loss  in  accuracy,  owing  to  the 
familiar  orthogonality  relations  (see  e.g.  Section  8.1  of  [Blaha,1977]). 

This  is  not  the  case  when  both  the  S.H.  and  P.M.  parameters  are  present, 
and  the  adopted  procedure  of  consecutive  adjustments  is  only  an  approxima¬ 
tion  to  a  simultaneous  solution.  However,  due  to  perhaps  insurmountable 
computer  requirements  accompanying  such  a  simultaneous  solution,  the  con¬ 
secutive  adjustment  algorithm  offers  a  viable  alternative  allowing  for  a 
great  reduction  in  the  number  of  parameters  solved  for  in  any  individual 
adjustment.  The  first,  global  (S.H.)  adjustment  leads  to  a  higher-order 
surface  (it  could  be  called  a  "S.H.  geoid").  This  surface  then  serves  as 
a  reference  surface  in  separate  local  (P.M.)  adjustments  in  the  areas  cover¬ 
ed  by  a  dense  net  of  observations  consisting  typically  of  satellite  alti¬ 
meter  data.  The  number  of  point  masses  and  their  location  including  their 
depth  below  the  earth's  surface  are  stipulated  beforehand.  The  model  equa¬ 
tions  for  both  AFj  and  AF2  are  based  on  the  expressions  for  the  disturbing 
potential  (T). 

After  a  L.S.  solution,  predictions  of  various  quantities  inves¬ 
tigated  in  this  chapter  may  be  of  interest.  They  are  given  by  any  of  the 
right-hand  sides  in  equations  (4.1),  but  with  the  quantity  "Lb"  absent.  The 
corresponding  variances  are  obtained  as  follows: 
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U1 


"1  S.H.  "1  ’ 


a2  =  A2  SP.M.  A2  • 
a2  =  aj  +  » 

where  the  variance-covariance  matrix  E  for  the  parameters  (S.H.  or  P.M. ) 
is  obtained  as  the  inverse  of  the  normal  equation  matrix  in  the  pertinent 
adjustment.  The  simple  addition  of  variances  yielding  the  total  variance 
in  the  last  equation  above  stems  from  the  approximation  allowing  for  an  in¬ 
dependent  treatment  of  the  two  kinds  of  parameters  by  the  adjustment  algo¬ 
rithm. 

The  expressions  for  the  fundamental  quantity,  T  ,  will  now  be 
recapitulated.  The  S.H.  formulation  appeared  in  (3.75)  essentially  as 

N 

T  =  (kM/r)  l  (a/r)n  AS(n)  ,  (4.4) 

n=2 

where  the  symbol  a  has  replaced  aQ  ,  the  symbol  k  (the  gravitational 
constant)  has  replaced  G  ,  and  AS(n)  has  been  given  in  (3.76).  The  P.M. 
formulation  at  point  i  can  be  adopted  from  (5.9)  of  [Blaha,1979],  giving 

^  =  I  (l/^jMWOj  ,  (4.5) 

j 

where 

£..  *  distance  between  the  1-th  (observation)  point 

and  the  j-th  point  mass, 

(kM),  = 


j-th  scaled  point  mass  magnitude. 


In  the  spherical  approximation,  R  represents  the  earth's  mean  radius 
(6371km)  and  G  is  the  average  value  of  gravity  (980  gals)  on  the  earth's 
surface;  the  point  masses  are  assumed  to  be  at  the  depth  d  below  the  sphere 
of  radius  R,  so  that 


-  (r2  +  Ri 


(4.6a) 


R1  =  R-d  , 

Fij  =  R  R1  cos  *ij  * 

cos  ip. .  =  sind).  sin<j>.  +  cosd> .  cos<j>.  cos(X. -X.)  , 

’  J  I  J  *  J  I  sJ 


(4.6b) 
(4.6c) 
(4 . 6d) 


where  41^  is  the  spherical  distance  between  the  points  i,j  with  the  coordinates 
(♦j.X^)  and  (4>j  *xj ) »  respectively. 

We  now  develop,  in  a  straightforward  fashion,  several  expressions 
where  the  differentiation  is  performed  with  respect  to  4^.,  X^  and  r^  (r^  is 
the  radial  distance  from  the  geocenter  to  point  i,  numerically  r^R): 


=  ( i  )  RRi  [cos<b.  s  i  n  4> .  -  sin<j>.  cos<J>.  cos(X.  -X.)], 

'  J  A  I  J  '  J  '  J 


etc.,  which  lead  to 


3Ti/3*1  -  l  [9(l/-eij)/9<l>i](kM)j 

j 

=  RRj  l  ( l/^ij-3)  tcosd)^  s1n<f>j  -  sin^  cos$j  cos(X.  -  X^)]  (kM)^  , 
J 


(4.7a) 
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(4.7b) 


3T1/3Xi  *  -RRj  COS^.  I  (1/£.J3)  cos^.  sin(X.  -  XjMkM)^  , 

J 

8Vari  =  -I  (R- cos^jXkN)^  ,  (4.7c) 

J 

32T./3<{.|  =  RRj  l  (1  /l.p  {(3RRj /l.p 
3 

«[cos<j> .  sin4> .  -  sin<j>.  cos*.  cos(X.  -  X  .)]2  -  cos^.  .}(kM) .  ,  (4.7d) 

1  j  1  j  '  j  ij  j 

32T./3X|  =  RRj  l  (1  /l.*)  (ORRj/^.p 
J 

* [cos<}>..  costj  sin(x.-x.)]2  -  cos^.  cos^.  cos(X.  -  XjOXkM) .  ,  (4.7e) 

dzT./Zrf  =  l  (1  /l.p  [3(R -  R1  cos ♦jjJVtjj2-  l](kM) ^  ,  (4.7f) 

J 

32Ti/a<bi 3Xi  =  -RRj  l  (1/^-jj3)  ((SRRj/^j2)  cos^  cos^  slnCx^  -  A^) 

3 

- [cosd> -  siruj).  -  simji.  cos<*>.  cos(X.  -  X.)J 

-  sirvjK  cos^.  s1n(X1  -  Xj»(kM)j  ,  (4.7g) 

a2Ti/3«j.i3ri  =  -Rj  l  (1/^j3)  [co$<f>.  sin^-sin^  cos^  cos(Xi -X.)] 

j  J 

«[3R(R-Rj  cos*fj)/*f j2  -  l](kH)j  ,  (4.7h) 

32T./3X.3r.  =  Rj  cos4>^  [  (l/£y3)  cos<(>j  sinUj-Xj) 

J 

■[3R(R-Rj  cos^j)//^2-  l] (kM) j  .  (4.7i) 
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The  mixed  second  derivatives  can  be  verified  upon  reversing  the 
order  of  differentiation.  Fur thermo re, the  Laplace  condition  should  be  ful¬ 
filled  by  construction.  This  is  verified  as  follows.  Upon  applying  (3. 73') 
at  point  i  we  have 


Txx  =  (1/R)  3Ti/3ri  +  (1/R2)  32V3«>i 


(4.8) 


=  l  (l/^1-j3){(3Rj/fij2)[cos<j)i  slityj-slrty.  cos<Jk  cos(Ai  -  A^)] 2  -  1  > ( kM ) ^ 
J 


and  from  T  and  Tzz  In  (3.73)  we  deduce 

Tyy  =  ?  (1/£ij3)  CC3Ri/^i/)  c°s2<frj  sin2(Ai-Xj)  -  l](kM)j 

3 


(4.9) 


Tzz  =  l  (1/^ij3)  [3(R'Ri  cos^j)2/^/-  l](kM)j.  , 
3 


(4.10) 


(4.10)  being  essentially  (4. 7 f ) .  From  (4.8)  -  (4.10),  straightforward  algebra 
yields 


T  +  T  +  T  =0 
xx  yy  zz 


(4.11) 


as  it  should.  For  reference  purposes,  from  (3,73)  we  also  have 


Txy  =  ‘3R1  £  (1/£i/)  cos<^j  sin(Ai  ‘ 

J 

■  [cos^  sin4> .  -  s in4> -  coscp .  cos(A.  -  A.)](kM) . 

*  J  *  J  *  J  J 


(4.12) 


which,  however,  will  not  be  used  in  the  balance  of  this  report. 
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In  the  following  sections,  five  observational  modes  will  be  in¬ 
vestigated:  geoid  undulations,  gravity  anomalies,  deflections  (north  and 
east)  of  the  vertical,  horizontal  gradients  (north  and  east)  of  gravity 
disturbance,  and  vertical  gradients  of  gravity  disturbance.  In  a  first 
step  for  each  mode,  the  S.H.  formulation  for  F(Xj,0)  and  AF^A^Xj  with  the 
row  matrix  Aj  given  explicitly  will  be  presented;  this  part  will  not  con¬ 
tain  the  familiar  spherical  approximation.  In  a  subsequent  step,  the  P.M. 
formulation  for  AF2  =  A2X2  with  A2  given  explicitly  will  be  considered  in 
spherical  approximation.  Also  included  will  be  the  practical  discussion 
pertaining  to  the  cut-off  distance,  l.e.,  to  the  size  of  a  spherical  cap 
centered  at  the  observation  point  beyond  which  the  P.M.  parameters  are 
Ignored  during  the  formation  of  an  observation  equation. 
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4.2  Geoid  Undulations 


For  the  most  part  this  model  has  been  already  treated  in  [Blaha, 
1977]  and  [Blaha, 1979] ,  referred  to  in  this  chapter  as  [Bl]  and  [B2] ,  re¬ 
spectively. 


S.H.  part.  Upon  the  customary  assumption  Wo  =  UQ,  stating  that 
the  potential  of  the  geoid  equals  the  potential  of  the  geocentric  reference 
ellipsoid  and  is  not  subject  to  adjustment,  we  have  from  Section  2.1  of 
[B2]: 


where 


with 


N  =  N(1  +  c)  , 


N 

N  =  rQ  l  (a/r1)"  AS(n)  , 


(4.13) 

(4.14) 


r 


o 


kM/Wo  , 


r1  =  a/[l  +  e2  sin2<f>/(l  -  e2)]'5  , 


r1  being  the  radial  distance  from  the  geocenter  to  the  ellipsoid  at  the 
point  where  N  is  being  sought,  a  and  e  being  the  semi -major  axis  and  the 
eccentricity  of  this  ellipsoid,  respectively,  and  $  being  the  geocentric 
latitude  of  the  ellipsoidal  point  in  question;  and  where 

c  =  Cj  sin2^  +  c2  ,  (4.15) 
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with 


cl  =  -3.0033  C2Q  -  0.005169  =  -0.00192  , 
c2  =  1.0011  C2Q  +  0.005169  *  0.00408. 

The  expression  (4.13)  corresponds  to  F(Xj,0)  upon  the  recognition  that  the 
initial  values  of  the  S.H.  parameters  enter  into  the  definition  of  AS(n). 

A  similar  statement  will  be  applicable  in  the  following  sections,  and  it 
need  not  be  repeated. 

If  the  geoid  undulations  are  part  of  a  satellite  altimetry  model, 
the  model  equation  becomes 

H  =  R  -  r  + (correction) 

as  described  in  Section  2.4  of  [Bl]  ,  where  H  is  the  distance  from  the  satel¬ 
lite  to  the  (S.H.)  geoid,  R  is  the  radial  distance  from  the  geocenter  to 
the  satellite,  and 

r  *  r'  +  N  ,  (4.16) 

(correction)  =  4.9m  sin2  2<p\.,  per  1000km  of  H  . 

The  parameters  X^  are  the  corrections  to  the  initial  values  of 
the  S.H.  potential  coefficients  Cno,  Cnn),  Snm  (referred  to  as  C's  and  S's), 
and  the  row  matrix  A1  Is  composed  of  the  partial  derivatives  of  N  with  re¬ 
spect  to  these  coefficients.  Most  often  the  value  c  in  (4.13)  could  be 
ignored  during  the  differentiation,  the  greatest  possible  relative  error 
thus  introduced  being  0.4%;  in  this  case  we  have 
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3N/9Cnm  =  r0^a/r'^n  cos  mA  Pnm^sin<t)^  ’ 

3N/3Snm  =  ro(a/r’)n  sin  mX  Pnm(sin$)  . 

For  more  accuracy  these  partial  derivatives  can  be  multiplied  by  (1  +  c). 

In  order  to  eliminate  the  remaining  error,  N(l. 0011  -  3.0033  sin2^)  can  be 
added  to  3N/3C2Q. 

If  used  as  a  part  of  the  satellite  altimetry  model,  the  above 
partial  derivatives  are  combined  with  those  pertaining  to  the  state  vector 
parameters  as  demonstrated  in  Sections  2.4  or  4.1  of  [Bl] .  In  the  latter 
section,  the  partial  derivatives  3N/3Cnm,  etc.,  are  replaced  by  the  equiva¬ 
lent  expressions  for  3r/3Cnm,  etc. 

P.M.  part.  As  in  Section  5.2  of  [B2] ,  the  part  AF,,  attributed 
to  the  point  masses  is 

Nn-  =  (1/S)  l  (1/^-jj)  (kM)j  ,  (4.17) 

J 

which  can  be  evaluated  with  the  aid  of  equations  (4.6).  The  row  matrix  A,, 
can  be  formed  directly  from  (4.17)  with  the  (kM).'s  absent  and  with  the  j's 

J 

determining  the  ordering  of  the  elements.  Since  all  the  investigated  ob¬ 
servational  modes  are  linear  insofar  as  the  P.M.  parameters  are  concerned, 
the  formation  of  Ag  is  thus  settled  once  and  for  all. 

An  important  decision  to  be  made  when  formulating  the  P.M.  ad¬ 
justment  model  in  practice  regards  the  cut-off  distance.  In  addressing 


this  task,  only  one  point  mass  (kM).  is  considered  and  the  position  of 

J 

the  observation  point  i  on  the  earth's  surface  (here  a  sphere  of  radius 
R)  is  varied  relative  to  this  point  mass.  According  to  [B2],  the  depth  of 
the  point  masses  is  taken  as  a  1.6 -multiple  of  their  horizontal  separa¬ 
tion  (as  projected  on  the  earth's  surface).  In  particular,  if  the  separa¬ 
tion  in  degrees  of  arc  is  s°  =  2°,  then 

s  *  2?2.4  km  , 
d  =  1.6s  =  350.0  km  , 

Rj  =  6021  km  . 

We  now  specialize  (4.17)  for  j  =  l,  giving 

Ni  *  (l/-ei:j)(kM)./G  ,  (4.17') 

where  the  one  point  mass  present,  ( kM) j  ,  is  chosen  to  be  a  10~6-part  of  the 
earth's  kM  (kM  =  3.986xl014  m3/sec2).  Next  the  position  of  the  observation 
point  i  is  varied  and  the  corresponding  value  computed.  This  varia¬ 
tion  can  proceed  by  s°  (with  intermediate  values  of  interest  included),  which 
makes  it  possible  to  see  at  how  many  s-intervals  away  from  the  observation 
point  the  P.M.  parameters  will  still  significantly  affect  the  modeled  value 
(and  vice  versa)  and  should  therefore  be  included  in  the  observation  equa¬ 
tion  or  in  the  prediction  formulas.  In  theory  one  could  include  all  of  the 
P.M.  parameters  in  every  observation  equation,  but  this  would  often  result 
in  prohibitive  computer  run-time  requirements.  The  specialization  illustrat¬ 
ed  in  (4.17')  is  typical  for  the  present  cut-off  analysis  and  will  be  used 


for  all  the  other  investigated  modes  as  well  (the  values  of  s,  d  and  will 
be  unchanged  from  the  above).  Therefore,  little  additional  explanation  will 
be  needed. 

The  procedure  just  described  leads  to  the  values  below  which  could 

serve  to  plot  a  curve  of  geoid  undulations  (N.)  generated  by  the  given  value 

(kM).  for  various  angles  ip.,  (in  degrees): 

J  1 0 


0°  ., 

. .  116  m 

6°  ... 

55.2  m 

2°  . 

, .  98,9  m 

8°  ... 

43.6  m 

4°  .. 

. .  73.1  m 

O 

o 

35.8  m 

The  undulation  curve  decreases  asymptotically.  At  the  distance  4s  (cor¬ 
responding  to  8°  above)  the  effect  of  the  P.M.  has  diminished  to  about  37.6% 
of  the  maximum  effect  exercised  by  a  P.M.  directly  underneath  the  observation 
point.  Compelled  by  practical  considerations  --  mainly  by  the  fact  that  the 
curve  decreases  rather  slowly  beyond  4s  --  we  can  adopt  this  distance  as 
the  cut-off  distance  with  regard  to  both  observations  and  predictions. 
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this  task,  only  one  point  mass  (kM).  is  considered  and  the  position  of 

J 

the  observation  point  i  on  the  earth's  surface  (here  a  sphere  of  radius 
R)  is  varied  relative  to  this  point  mass.  According  to  [B2] ,  the  depth  of 
the  point  masses  is  taken  as  a  1.6 -multiple  of  their  horizontal  separa¬ 
tion  (as  projected  on  the  earth's  surface).  In  particular,  if  the  separa¬ 
tion  in  degrees  of  arc  is  s°  =  2°,  then 

$  a  222.4  km  , 
d  =  1.6s  =*  350.0  km  , 

Rj  =  6021  km  . 

We  now  specialize  (4.17)  for  j  =  l,  giving 

N.  =  (l/^ij)(kM)J/G  ,  (4.17*) 

where  the  one  point  mass  present,  (kM).  ,  is  chosen  to  be  a  10~6-part  of  the 

J 

earth's  kM  (kM  =  3.986xl014  m3/sec2).  Next  the  position  of  the  observation 
point  i  is  varied  and  the  corresponding  value  computed.  This  varia¬ 
tion  can  proceed  by  s°  (with  intermediate  values  of  interest  included),  which 
makes  it  possible  to  see  at  how  many  s-intervals  away  from  the  observation 
point  the  P.M.  parameters  will  still  significantly  affect  the  modeled  value 
(and  vice  versa)  and  should  therefore  be  Included  in  the  observation  equa¬ 
tion  or  in  the  prediction  formulas.  In  theory  one  could  include  all  of  the 
P.M.  parameters  in  every  observation  equation,  but  this  would  often  result 
in  prohibitive  computer  run-time  requirements.  The  specialization  illustrat¬ 
ed  in  (4.17')  is  typical  for  the  present  cut-off  analysis  and  will  be  used 
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4.3  Gravity  Anomalies 


As  in  the  previous  section,  much  of  the  development  in  [Bl]  and 
[B2]  will  now  be  briefly  recapitulated  and  complemented  by  the  cut-off  con¬ 
siderations. 

S.H.  part.  A  satisfactory  yet  simple  model  has  been  analyzed 
in  Chapter  5  of  [Bl] .  From  (5.24a)  of  this  reference  we  write 

N 

Ag  *  C  l  (n-l)(a/r)n  AS(n)  ,  (4.18) 

n=2 

where 

C  =  kM/r2 

with  r  obtained  as  in  (4.16). 

According  to  equations  (5.27)  of  [Bl]  ,  Aj  is  composed  of  the 

elements 

9Ag/3Cno  =  C(n-l)(a/r)n  Pn(sin^)  , 

3Ag/8Cnm  =  C(n-l)(a/r)n  cos  mX  Pnm(sin$)  , 

3Ag/3Snn)  =  C(n-l) (a/r)n  sin  mX  Pnm(sin$). 

The  gravity  anomaly  model  as  well  as  the  above  partial  derivatives  have  been 
refined  in  Chapter  4  of  [B2] .  Exceptionally,  such  a  model  could  be  used  for 
maximum  accuracy,  but  there  is  no  need  to  rewrite  all  the  pertinent  formulas. 
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P.M.  part.  According  to  (5.10)  of  [B2]  ,  we  have 

Agi  =  (l/R)  l  (l/^ij)[(R2-Fij)/£ij2-  2](kM)j  ,  (4.19) 

3 

which,  as  in  the  case  of  geoid  undulations,  can  be  evaluated  by  means  of 
equations  (4.6). 

The  results  pertaining  to  the  cut-off  distance  are  presented  as 

follows: 

0°  ...  290  mgal  6°  ...  21.2  mgal 

2°  ...  172  mgal  10°  ...  1.0  mgal 

4°  ...  62.0  mgal  20°  ...  -3.0  mgal 

Clearly,  the  cut-off  distance  could  be  chosen  as  2s  (here  4°)  where  the  ef¬ 
fect  of  the  P.M.  has  diminished  to  about  21.4%  of  the  maximum  effect  exercis¬ 
ed  by  a  P.M.  underneath  the  observation  point.  The  3s-cap  (here  6°)  would  be 
more  than  sufficient,  the  corresponding  effect  shrinking  to  a  mere  7.3%. 
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4.4  Deflections  of  the  Vertical 


A  component  of  the  deflection  of  the  vertical  (e)  in  any  azimuth 
is  given  (see  e.g.  [Heiskanen  and  Mori tz, 1967] ,  page  112)  as 

e  =  -dN/ds  . 

In  the  north-south  direction  we  have 
£  =  -dN/ds j  , 

and  in  the  east-west  direction, 
n  =  -dN/ds2  , 

where  dSj  and  ds2  are  positive  with  increasing  latitude  and  longitude,  re¬ 
spectively,  and  dN  is  the  change  in  geoid  undulation  along  the  appropriate 
directions.  From  the  geometry  of  a  spheroid  It  follows  that 

dSj  =  p  d<j>  =  r'  d?  (1  -  2e2  cos  2$)  , 

ds2  =  v  cos<p  dA  =  r'  cos?  dA  , 

where  p  and  v  are  the  radii  of  curvature  in  the  meridian  and  in  the  prime 
vertical,  respectively.  It  thus  follows  that 

C  =  - ( 1/r ' ) ( 1  +  2e2  cos  2?)  9N/9?  ,  (4.20) 

n  =  -  [l/(r'  cos?)]  3N/9A  ,  (4.21) 


the  last  expression  not  to  be  used  for  the  poles.  In  the  spherical  ap¬ 
proximation  r1  would  be  replaced  by  R  and  e  would  be  replaced  by 
zero,  and  the  formula  (2-204)  of  the  above  reference  would  be  recovered. 


S.H.  part.  From  (4.14)  it  follows  that 

N 

(l/r')9N/3*=  (r  /r')  l  (a/r '  )n  AS(n)  , 

n=2 

where  AS(n)  has  been  defined  in  (3.78c).  Formulas  (4.13)  and  (4.20)  thus 
result  in 

N 

C  -  -(r  /r')(l  +  c)(l  +  2e2  cos  2^)  £  (a/r* )"  AS(n)  ;  (4.22) 

0  n=2 


in  spherical  approximation  this  would  reduce  to 

N 

C  *  -  I  AS(n)  . 
n=2 


Similarly,  (4.13),  (4.14)  and  (4.21)  yield 


n  =  -(l/cos<|>)(r  /r')(l  +  c)  l  (a/r')n  AS'(n)  ,  (4.23) 

0  n=2 


where  S ' ( n)  has  been  defined  in  (3.78a).  In  spherical  approximation  one 
would  write 

N 

n  -  -(l/cos<}>)  l  AS'(n)  . 
n=2 


When  forming  the  elements  of  for  it  follows  from  (4.22)  that 


3C/3Cno  =  -tt'(a/r')n  dPn(s1ii*)/d*  , 

3C/3cnm  =  -tt'(a/r')n  cos  mX  dPnm(sin$)/d$  , 

3^3Snm  =  sin  dPnm(sin$)/d^  , 

where 

t  =  (r0/r')(l  +  c)  , 

t'  =  (l  +  2e2  cos  2$)  . 

Similarly,  from  (4.23)  we  obtain 

an/acno  *  o, 

3n/3Cnm  =  (t/cos*)(a/r-)n  m  sin  mX  Pnm(sin$)  , 

3n/3snm  s  -'(t/cos<j>)(a/r’  )n  m  cos  mX  Pnm(sin$)  . 

In  spherical  approximation  all  of  t,  t'  and  (a/r')n  would  be  replaced  by 
unity. 

P.M.  part.  From  (4.20)  and  (4.21)  we  write  in  spherical  approxima 

tion: 

^  =  - [1/(GR)]  3T i / 3<J>i  , 
ni  =  -  [l/(GR  cos^)]  3T./3X.  , 


where  we  have  used 


Thus,  upon  utilizing  (4.7a)  and  (4.7b)  one  obtains 


5i  =  -oyG)  2  (^ij3) 

j 

-  [cos<p.  sin<j)j  -  s i n<p -  cos^.  cos(Xi  -  X^)]  ( kM) ^  ,  (4.24) 


n.  -  (Rj/G)  l  (!/£,*/)  cos<}>j  sin(X.  -  Xj)(kM)j 
0 


(4.25) 


With  regard  to  the  cut-off  decision,  either  (4.24)  or  (4.25)  is 
applied  in  conjunction  with  the  given  point  mass  as  usual.  It  is  useful  to 
carry  out  this  analysis  in  a  "canonical  form"  where  the  P.M.  latitude  is 
chosen  to  be  zero  and  where  the  position  of  the  observation  point  varies 
either  in  the  northerly  direction  (in  conjunction  with  ^ )  or  in  the  easterly 
direction  (in  conjunction  with  r^).  The  same  results  are  obtained  for  either 
quantity  as  will  become  clear  by  inspection.  They  are  listed  as 


0° 

...  0.0" 

6°  ... 

13.2" 

1° 

...  18.0" 

o 

OO 

CO 

■S-J 

2° 

...  25.4" 

10°  ... 

6.0" 

2.3° 

...  25.7" 

12°  ... 

4.4" 

4° 

...  20.5" 

20°  ... 

1.7" 

At  a  point  above  the  P.M.  the  thus  generated  deflection  is  zero  (a  result 
highly  plausible).  The  deflection  values  rise  sharply  until  slightly  more 
than  Is  (here  2°)  and  then  fall  off  quite  rapidly;  the  curve  joining  such 
values  approaches  zero  asymptotically.  The  P.M.  has  the  most  significant 
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effect  for  distances  between  0  and  3s  (here  6°),  At  the  distance  5s  (here 
10°),  which  can  be  adopted  as  a  cut-off  distance,  this  effect  has  diminish¬ 
ed  to  only  about  23%  of  the  maximum  effect  around  Is.  The  cut-off  distance 
of  4s  corresponding  to  about  34%  of  the  maximum  effect  could  also  be  adopt- 


4.5  Horizontal  Gradients  of  Gravity  Disturbance 


The  gradients  (horizontal  and  vertical)  of  gravity  can  easily  be 
reduced  to  the  gradients  of  gravity  disturbance  as  outlined  in  the  previous 
chapter.  Thus  W  ,  etc.,  can  be  transformed  to  T  ,  etc.,  where 

Z  X  Z  X 

T  =  W  -  U  ,  etc.  Another  possibility  is  to  obtain  T  ,  etc.,  directly 

ZX  t  a  ZX  ZX 

from  a  gradiometer  as  indicated  in  [Moritz, 1975].  Although  working  with 
TZ"X",  etc.,  would  be  more  rigorous,  the  errors  introduced  by  replacing  the 
triad  i"  by  the  triad  i  are  deemed  negligible  for  the  present  purpose.  In 
any  event,  replacing  the  formulas  (3.73)  by  (3.72),  if  needed,  is  an  easy 
matter.  When  dealing  with  the  gradients  of  gravity  disturbance  in  the  form 
Tzx,  etc.,  we  are  in  fact  using  a  "spherical  approximation"  but  only  insofar 
as  the  formulas  for  the  gradients  are  concerned;  the  value  of  r  in  the 
S.H.  part  is  not  replaced  by  R  or  by  another  constant  but  is  computed  as 
in  (4.16),  with  a  height  above  the  earth's  surface  eventually  added  in  case 
of  aerial  or  satellite  observations.  In  the  P.M.  part,  r  will  be  replaced 
by  R  as  usual . 

S.H.  part.  From  the  3rd  and  5th  relations  in  (3.73)  in  combina¬ 
tion  with  the  appropriate  expressions  of  (3,77),  we  deduce  for  the  hori¬ 
zontal  gradients  of  gravity  disturbance  in  the  north-  and  east-  directions, 
respectively: 

N 

T7¥  =  -C  l  (n+2)(a/r)n  AS(n)  ,  (4.26) 

**  n=2 

N 

T  =  -(C/cos 4>)  l  (n+2)(a/r)n  AS'(n)  ,  (4.27) 

y  n=2 
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where 


I  C  =  kM/r3  , 

and  where  AS(n)  and  AS'(n)  have  been  defined  in  equations  (3.78). 

The  row  matrix  pertaining  to  the  "north"  gradient  is  seen 
from  (4.26)  to  be  composed  of  the  elements 

3Tzx/3Cno  =  -C(n+2)(a/r)n  dPn(sin$)/d$  , 

3Tzx/3Cnm  =  "C(n+2)(a/r)n  cos  mA  dPnm(sin?)/<4  , 

3Tzx/3Snm  =  -C(n+2)(a/r)n  sin  mA  dPnm(sin^)/d^  . 

Similarly,  with  regard  to  the  "east"  gradient  we  have  from  (4.27): 

3T  u/3Crt  -  0  , 
zy  no 

3Tzy/3Cnm  =  (C/cos^) (n+2)(a/r)n  m  sin  mA  Pnm(sin^)  , 

3Tzy/3Snm  =  "(C/cos^)(n+2)(a/r)n  m  cos  mX  Pnm(sin^)  . 

P.M.  part.  Upon  replacing  r  by  R  in  the  3rd  and  5th  relation 
in  (3.73)  we  write  for  the  observation  point  i  in  the  P.M.  model: 

■  Tzx  =  (l/R)32T1/3<J>i3ri  -  (l/R2)3Ti/3<bi  > 

: 

Tzy  =  [l/(R  COS*i)]32T1/3Xi3r1  -  [l/R2  cos<}>. )]  3T./3Xi  . 
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With  the  use  of  the  appropriate  expressions  from  (4.7)  this  becomes 


Tzx  =  "3Ri  l  (i/^i j5)  (R-Rj  cos^) 

J 

«  [cos^.  sin4> -  -  sin<j>.  cost)).  cos(X.  -  A.)]  (kM) .  ,  (4.28) 

•  J  *  J  ’  u  J 


Tzy  =  3R1  E  (1/^i  ?)  (R  -  Rj  cos^j)  cos<jK  sin(X.  -  X^)(kM)^ 
j 


(4.29) 


It  can  be  noticed  that  the  parts  beyond  (R-R1  cosijj.,)  in  these  equations 

A  I  J 

agree  with  the  corresponding  parts  in  (4.24),  (4.25). 

Similar  to  the  analysis  carried  out  for  the  deflections  of  the 

vertical,  the  results  obtained  for  T  and  T  in  the  "canonical  form"  are 

zx  zy 

again  indentical.  Listed  below  are  a  few  values  in  E  (Eotvos,  0.1  mgal/km 
or  10"  sec"  )  obtained  by  this  procedure: 


0°  ...  0.00  E 
1°  ...  6.68  E 
1.8°  ...  7.75  E 
2°  ...  7.54  E 
4°  ...  3.44  E 


6°  ...  1.33  E 
8°  ...  0.58  E 
10°  ...  0.29  E 
12°  ...  0.17  E 
20°  ...  0.04  E 


A  curve  constructed  from  these  values  would  be  steeper  than  its  counterpart 
for  the  deflections  of  the  vertical.  Again,  the  gradient  generated  above 
the  P.M.  is  zero.  The  P.M.  has  the  most  significant  effect  between  0  and  2s 
(here  4°).  The  cut-off  decision  could  be  made  in  favor  of  3s  (here  6°)  where 
this  effect  has  diminished  to  about  17%  of  the  maximum  effect  around  Is. 


4.6  Vertical  Gradients  of  Gravity  Disturbance 


S.H.  part.  From  the  last  relation  in  (3.73)  and  from  the  last 
expression  of  (3.77)  we  obtain  the  vertical  gradient  of  gravity  disturbance 
as 

N 

T7  =  C  l  (n+l)(n+2)(a/r)n  AS(n)  ,  (4.30) 

n=2 

where  AS(n)  has  been  defined  in  (3.76), 

Accordingly,  contains  the  elements 

9Tzz/9Cno  =  C(n+l)(n+2)(a/r)n  Pn(sin$)  , 

9Tzz/9Cnm  =  C(n+l)(n+2)(a/r)n  cos  mX  Pnm(sin4>)  , 

9Tzz/9Snm  s  C(n+l)(n+2)(a/r)n  sin  mX  Pnm(sinc|>)  . 


P.M.  part.  In  spherical  approximation,  the  last  relation  of 
(3.73)  is  written  for  the  observation  point  i  in  the  P.M.  model  as 


T 


zz 


a^./ar? 


9 


which,  according  to  (4.7f),  is 

Tzz  =  2  d/^ij3)  [3(R-RX  cos^j)2/^./-  l](kM)..  .  (4.31) 

j 


-87- 


The  procedure  involving  the  "canonical  form"  and  one  point  mass 

results  in 

0°  ...  18.59  E  6°  ...  -0.19  E 

2°  ...  6.97  E  8°  ...  -0.21  E 

4°  ...  0.67  E  10°  ...  -0.15  E 

The  vertical  gradient  curve  is  thus  seen  to  be  steeper  than  its  counterpart 
for  the  horizontal  gradients.  The  maximum  effect  corresponds  now  to  ^ =  0 . 

*  J 

The  cut-off  distance  of  2s  (here  4°)  is  more  than  adequate,  the  effect  of  the 
P.M.  at  that  distance  being  about  3,6%  of  the  maximum  effect  (at  Is  this  ef¬ 
fect  would  correspond  to  37.5%). 


5.  SPECIFIC  CONSIDERATIONS  RELATED  TO 
SEA  SURFACE  TOPOGRAPHY 


5.1  Outline 


The  sea  surface  topography  represents  a  dynamic  surface  whose 
slight  changes  in  shape  are  imputable  to  astronomical  forces  in  combina¬ 
tion  with  other  effects  (oscillations  in  water  basins,  shallow  water  ef¬ 
fects,  meteorological  effects,  etc.).  The  astronomical  forces  are  at 
the  basis  of  tidal  changes  in  the  sea  surface.  The  attraction  of  the 
moon  represents  the  most  important  aspect  of  these  forces.  By  comparison, 
the  effect  of  the  sun  amounts  to  only  about  46%  of  the  effect  of  the  moon, 
both  effects  having  similar  characteristics. 

If  the  moon's  orbit  were  in  the  plane  of  the  earth's  equator, 
there  would  be  two  high  tides  and  two  low  tides  in  a  lunar  day  (approxi¬ 
mately  24  hours  and  50  minutes,  or  1.035  solar  day),  and  the  tides  would 
be  semi-diurnal.  If  the  plane  of  the  ecliptic  coincided  with  the  plane 
of  the  earth's  equator,  the  sun  would  produce  similar  effects  except  for 
the  magnitude  and  the  period  of  the  tidal  cycle  (a  solar  day  has  24  hours). 
The  varying  declinations  of  the  moon  and  the  sun  produce  additional  ef¬ 
fects  that  are  essentially  diurnal.  Irregularities  of  the  moon's  orbit 
around  the  earth  and  in  the  earth's  orbit  around  the  sun  (varying  dis¬ 
tance  and  varying  velocity)  require  additional  considerations  as  explain¬ 
ed  e.g.  in  Chapter  2  of  [Macmillan,  1966]. 
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Further  complications  are  due  to  the  departure  of  the  earth's 
structure  from  an  idealized  model  which  would  give  rise  to  the  "equili¬ 
brium  tide".  Assumptions  associated  with  the  equilibrium  tide  include  a 


completely  rigid  earth  covered  over  its  entire  surface  by  deep  water, 
absence  of  friction  in  the  water  envelope,  absence  of  meteorological  and 
other  disturbances,  etc. 

In  reality,  tidal  changes  include  two  important  effects  expres¬ 
sed  by  Love's  numbers;  these  effects  are  the  ocean  floor  deformation  and 
an  additional  tide  induced  by  such  a  deformation.  If  included  in  the 
tidal  model,  they  would  result  in  modifying  the  amplitude  of  the  equili¬ 
brium  tide  to  a  limited  extent.  However,  this  consideration  is  not  al¬ 
ways  of  paramount  importance  in  applications  of  satellite  altimetry.  In 
particular,  when  one  proceeds  to  adjust  satellite  altimetry  in  a  localized 
area  this  amplitude  may  also  be  subject  to  adjustment  and  the  inclusion 
of  a  factor  due  to  the  above  two  effects  would  simply  lead  to  the  replac¬ 
ing  of  one  unknown  parameter  by  another  (the  altimeter  does  not  offer 
means  for  separate  solution  of  Love's  numbers).  Accordingly,  only  the 
equilibrium  tide  will  be  considered  in  this  discussion. 

In  the  satellite  altimetry  mathematical  model  the  ocean  surface 
is  considered  as  equipotential .  At  a  given  moment  the  equilibrium  model 
results  in  an  equipotential  surface  (under  the  assumption  of  an  idealized 
spherical  earth  it  is  a  spheroid).  Time  averaging  applied  to  this  surface 
also  results  in  an  equipotential  surface  (a  different  spheroid).  Accord¬ 
ingly,  if  the  altimeter  could,  over  a  sufficiently  long  time  span,  sense 
the  latter  (average  equipotential)  surface,  the  two  models  would  indeed 


be  compatible.  However,  this  Is  not  the  case  because  of  several  effects 
not  accounted  for  Including  water  friction,  continent  boundaries,  currents, 
meteorological  factors,  etc.  Consequently,  the  surface  sensed  over  a 
given  time  span  exhibits  slopes  in  different  areas  with  respect  to  anequi- 
potential  surface.  If  an  adjustment  does  not  allow  for  the  inclusion  of 
such  slopes,  the  parameters  describing  the  equlpotential  surface  may  be 
slightly  deformed  so  as  to  accommodate,  in  a  least-squares  sense,  the  mea¬ 
sured  surface.  It  is  thus  desirable  to  provide  for  additional  parameters 
describing  the  sea  surface  slopes,  at  least  in  the  areas  where  the  slopes 
are  suspected  or  known  to  exist.  If  the  time  history  of  such  slopes  is  of 
interest  (assuming  that  the  slopes  are  changing)  the  adjustment  may  be  re¬ 
peated  at  suitable  time  intervals  provided  sufficient  altimeter  data  are 
available.  A  simple  model  will  be  developed  shortly  that  allows  for  the 
inclusion  of  new  parameters  in  an  adjustment,  accommodating  the  slopes 
(including  a  vertical  separation)  between  the  measured  surface  and  the 
equi potential  surface  in  predetermined  areas. 

However,  prior  to  this  task  we  shall  direct  some  attention  to 
the  problem  of  a  reference  surface  for  heights,  which  we  shall  call  the 
fixed,  or  static,  geoid.  In  [Bomford,  1975],  page  247,  it  is  suggested 
that  "the  geoid  (and  the  datum  of  height)  is  an  equipotentlal  surface  of 
the  earth's  attraction  and  centrifugal  force  only,  not  including  the  sun 
and  moon...".  It  is  also  stated  that  "for  its  datum  surface  for  height, 
every  country  or  group  of  countries  select  mean  sea-level  at  a  defined 
tide  gauge...".  However,  even  if  the  best  possible  mean  (in  time)  sea 
level  were  selected  at  given  points,  it  could  not  represent  globally  a 


consistent  datum  for  heights,  i.e.,  the  static  geoid.  This  can  be  visu¬ 
alized  as  follows.  Suppose  that  the  static  geoid  is  a  perfect  sphere 
with  a  known  center  and  of  the  same  volume  as  a  mean  sea  surface  cover¬ 
ing  the  whole  globe,  assuming  the  equilibrium  model.  Due  to  the  moon's 
and  sun’s  attractive  forces,  the  mean  sea  surface  is  a  spheroid.  If  one 
wishes  to  define  the  static  geoid  (a  sphere  in  this  context)  passing 
through  a  given  tide  gauge  on  the  equator,  one  is  faced  with  the  problem 
that  this  sphere  is  larger  than  that  representing  the  static  geoid  pas¬ 
sing  through  a  given  tide  gauge  at  a  higher  latitude.  Thus,  in  order  to 
define  a  unique  global  datum  for  heights,  a  suitable  correction  should 
be  applied  to  the  mean  sea  level  as  a  function  of  latitude.  This  problem 
is  discussed  in  the  following  section  where  the  moon's  effect  is  consider¬ 
ed  first  and  the  moon-earth  distance  is  assumed  constant,  equal  to  the 
mean  distance.  The  sun's  effect  is  incorporated  along  similar  lines. 
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5.2  Considerations  Related  to  Static  Geoid 

In  this  section  the  equipotential  surface  without  the  influence 
of  astronomical  forces  is  assumed  to  be  a  sphere.  This  would  be  equiva¬ 
lent  to  considering  the  earth  a  point  mass  or  a  homogeneous  nonrotating 
sphere,  or  simply  to  retaining  only  the  leading  term  in  the  expansion  of 
the  earth's  potential  in  spherical  harmonics,  etc.  In  determining  the 
changes  in  the  potential  due  tothemoon  and/or  sun  such  approximation  is 
acceptable.  The  static  geoid  is  thus  considered  to  be  a  sphere  and  the 
main  task  consists  in  determining  its  radius  (R)  or,  equivalently,  its  re¬ 
lation  to  the  mean  sea  surface.  The  assumption  of  the  equilibrium  tide 
is  used  throughout. 

In  considering  the  moon's  effect  In  this  model,  the  total  poten¬ 
tial  at  a  point  (P)  on  the  static  geoid  is  determined  as  the  sum  of  the 
potential  generated  by  the  earth's  mass  concentrated  at  its  center,  of  the 
potential  due  to  the  moon's  mass  (concentrated  at  the  moon's  center)  and 
of  the  potential  generated  by  the  monthly  rotation  of  P  around  the  bary- 
center  of  the  earth-moon  system.  A  new  equipotential  surface  passing 
through  P  Is  computed  (see  [Bomford,  1975],  pp.  271,  272)  through  its 
"tidal  undulation",  N,  with  respect  to  the  sphere  of  radius  R.  In  this 
process,  the  volume  enclosed  by  the  two  surfaces  is  constrained  to  be  the 
same,  i.e.,  the  mean  value  of  N  over  the  sphere  is  constrained  to  zero. 
This  is  plausible  since  the  volume  of  the  water  envelope  as  well  as  the 
volume  of  the  solid  earth  beneath  It  do  not  change  —  no  mass  is  added  or 
subtracted.  Upon  considering  the  moon's  orbit  as  circular  and  adopting 


the  average  earth-moon  distance  for  this  purpose,  the  tidal  undulation  N 
for  any  time  and  place  is 


N  *  k(cos  2z  +  1/3)  ,  (5.1) 

where  z  is  the  angle  between  the  directions  toward  P  and  toward  the  moon 
measured  at  the  earth's  center;  it  is  approximately  equal  to  the  zenith 
distance  of  the  moon  at  P  (the  angle  under  which  the  earth's  radius  is 
seen  from  the  center  of  the  moon  is  neglected).  The  angle  z  is  computed 
as  a  function  of  geographic  location  (of  P)  and  of  time.  The  above  ap¬ 
proximation  of  circular  orbit  yields 

k  «  kj  =  0.267  m .  (5.2a) 

In  case  of  the  sun  (5.1)  would  hold  again,  but  with 

k  =  k2  =  0.123m,  (5.2b) 

which  indicates  that  the  tidal  amplitude  associated  with  the  sun  is  about 

46%  of  that  associated  with  the  moon.  If  an  overbar  denotes  the  mean 
over  a  unit  sphere,  it  is  confirmed  that  N  =  Q  because  at  a  given  moment 
cos  2z * -1/3. 

It  will  be  shown  shortly  that  the  surface  determined  by  (5.1) 

Is  a  prolate  spheroid  (ellipsoid  of  revolution)  with  its  long  axis  point¬ 
ing  toward  the  moon.  This  situation  is  schematically  illustrated  in 
Figure  5.1  where  the  "z"  corresponding  to  P j  ,  P2  and  P'  could  be  denoted 
as  Zj  ,  z2  and  z'  ,  respectively;  in  particular,  Zj  =  0,  z2  =  90°  and 


z'  =  54.7°  .  (5.3) 

(There  are  four  such  values  of  z*  corresponding  to  N=0,  symmetric  about 
the  axes  x,y.)  Associated  with  (5.2a),  we  have 

Nj  =  0.356m,  N2  =  -0.178m,  (5.4) 

while  for  (5.2b)  these  undulations  would  he  0.164m  and  -0.082m,  respec¬ 
tively;  if  both  the  moon  and  sun  were  in  line  with  the  earth  (i.e. ,  at 
conjunction  or  opposition)  the  corresponding  values  would  be  added  alge¬ 
braically  and  the  tidal  undulations  would  reach  0.52m  and  -0.26m,  re¬ 
spectively.  If,  in  this  case,  the  x  axis  in  Figure  5.1  depicted  the 
earth's  equator,  the  tidal  range  on  the  equator  would  reach  about  0.78m 
(i.e.,  0.52m  +0.26m)  ;  this  order  of  magnitude  is  in  agreement  with  the 
realistic  tidal  range  in  open  ocean. 

In  order  to  show  that  the  curve  representing  "spheroid"  in  Fig¬ 
ure  5.1  is  indeed  an  ellipse  —  or  very  nearly  so  —  we  adopt  from  the 
figure  and  from  (5.1): 

a  =  R  +  (4/3) k  ,  b  *  R  -  (2/3)k  ; 

upon  attributing  to  the  point  Q  the  coordinates  x=  (R  +  N)  cos  z  and 
y=(R  +  N)  sin  z  with  N  given  in  (5.1),  after  some  manipulation  we  obtain 

x2/a2  +  y2/b2  *  1  +  t  , 

where 

i  •  0.25  x  10-”  , 


-95- 


1 

i 

j 


Figure  5.1 

Schematic  representation  of  the  static  geoid  and  equilibrium  tide 


which  indicates  that  to  any  reasonable  approximation  the  curve  in  question 
is  an  ellipse. 

In  Order  to  acquire  a  good  idea  about  the  separation  between  the 
two  surfaces  depicted  in  Figure  5.1,  we  compute  an  "RMS  separation"  as  fol¬ 
lows: 

N2  =  k2  [cos22z  +  (2/3)cos  2z  +  1/9]  , 

-  _!  it  2ir 

cos22z  =  (47t)“  /  /  cos22z  sin  z  dz  du  . 

z=0  u=0 
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The  last  expression  eventually  yields  14/30;  with  cos  2z  given  earlier, 
we  find 


(N2)*5  =  0.596  k  . 

Repeating  also  a  previous  result,  we  have  for  the  effect  of  the  moon: 

N*0  ,  (5.5a) 

(N2)*5  =  0.16  m.  (5.5b) 


If  the  effect  of  both  the  moon  and  sun  at  conjunction  or  opposition  are 
considered,  the  value  in  (5.5b)  is  replaced  by  0.23m. 


The  main  point  of  the  present  discussion  is  the  consideration  of 
mean  values  over  a  given  time  interval.  The  (approximate)  zenith  distance 
z  is  changing  with  time  due  to  the  earth's  rotation  and  the  moon's  orbit, 
and  it  can  be  expressed  from  the  following  standand  formula  of  spherical 
astronomy  relating  the  horizon  and  hour  angle  systems  (see  e.g.  [Mueller, 
1969],  p.  38): 


*\ 


sin  z  cos  A  =  sinS  cos<j>  -  cos6  cos  h  sin4> 
sin  z  sin  A  =  -cos6  sin  h  > 

cos  z  =  sinfi  sin<J>  +  cosfi  cos  h  cos<p  , 


(5.6) 


where  A  is  the  azimuth  of  the  celestial  bo4y  (here  the  moon),  h  and  6  are 
Its  hour  angle  and  declination,  respectively,  and  <|>  is  the  (geocentric) 
latitude  of  the  observer.  From  (5.6)  it  readily  follows  that 
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cos  2z  =  -cos  2<J>  +  sin  2$  sin  26  cos  h  +  cos  2<p  cos26 
+cos  2<p  cosz6  cos2h  -  cos26  sin2h  . 


(5.7) 


With  regard  to  h,  great  simplifications  will  take  place  if  the  averaging 
is  carried  out  over  an  integer  number  of  lunar  days  (i.e. ,  over  an  integer 
number  of  2tt). 

The  time-averaging  process  for  6  is  facilitated  by  assuming  that 
the  plane  of  the  moon’s  orbit  coincides  with  the  ecliptic  which  makes  an 
e  =  23.5°  angle  with  the  equator.  In  reality,  the  plane  of  the  moon's  or¬ 
bit  is  inclined  by  about  5°  with  respect  to  the  ecliptic  and  thus  inter¬ 
sects  the  plane  of  the  equator  at  angles  ranging  from  a  minimum  of  18.5° 
to  a  maximum  of  28.5°.  However,  it  will  be  shown  that  this  approximation 
results  in  an  error  in  tidal  undulation  of  less  than  2cm  in  the  worst  case. 
Consequently,  the  formula  giving  the  declination  of  the  moon  can  be  tran¬ 
scribed  in  our  context  from  the  formula  below,  giving  the  declination  of  the 
sun.  In  particular,  the  relationship  between  the  right  ascension  and  eclip¬ 
tic  systems  yields  for  6=0  (the  ecliptic  latitude  of  the  sun): 

sinfi  *  sinX  sine  ,  (5.8) 

where  X  is  the  sun's  ecliptic  longitude  measured  counterclockwise  from  the 
vernal  equinox  (point  of  intersection  between  the  plane  of  the  ecliptic  and 
that  of  the  equator).  This  follows  from  equation  (3.11)  of  [Mueller,  1969]. 
We  only  have  to  remember  that  when  adopting  this  equation  for  the  moon,  X 
covers  the  interval  2tt  in  a  sideral  month  (27.32  days).  When  averaging,  it 
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will  be  advantageous  to  integrate,  with  respect  to  X,  over  an  integer  num¬ 
ber  of  7T  corresponding  to  an  integer  number  of  13.66  days.  If  the  varia¬ 
tions  in  e  were  significant,  the  integration  could  probably  be  extended 
only  over  13.66  days  which  would  correspond  to  tt  in  X  and  to  about  13x2ir 
in  h  (more  precisely,  to  13.2x27t  in  h).  However,  since  e  has  been  fixed 
in  our  model,  the  integration  can  be  extended  over  any  number  of  13.66  days 
which  will,  of  course,  include  a  very  large  number  of  1.035  days  (the  larger 
this  number  the  lesser  the  averaging  error  if  this  number  is  not  a  perfect 
integer).  This  device  also  allows  for  a  simple  mathematical  combination 
of  the  moon's  and  sun's  effects  not  only  because  the  "e"  are  identical  but 
also  because  in  the  case  of  the  sun  the  integration  over  merely  ir  in  A  in¬ 
volves  the  time  interval  of  half  a  year  (which  is  an  order  of  magnitude 
greater  than  13.66  days).  The  useful  outcome  of  this  argument  is  that  our 
simplified  —  although  sufficiently  accurate  --  model  allows  for  a  simulta¬ 
neous  time-averaging  over  an  integer  number  of  2tt  in  h  and  an  integer  number 
of  tt  in  X  for  both  the  moon's  and  sun's  effects. 


When  averaging  (5.7)  over  the  above  intervals,  the  second  term  on 
the  right-hand  side  vanishes  because  of  cos  h.  We  then  substitute  (5.8)  in 
(5.7)  and  deduce: 


Xj+nir 

(mr)"1 

f  (1  -  sin2e 

s1n2X)dX  =  1  -  4sin2e 

X1 

(2nrrr)“1 

hn+2mTT 

/*  fcoszh 

hJ  lsi"2h. 

j  dh  =  i*  , 

» 


9 


yielding 


<cos  2z>  *  -cos  24>  +  (1  -  15Sin2e)  [(3/2)cos  24>  -  %  ]  , 

where  <> indicates  a  time  average.  After  little  manipulation,  from  (5.1) 
we  obtain 

<N>  =  %k(cos  24)-  1/3) [l  -  (3/2)sin2e]  ,  (5.9) 

where  k  may  represent  k^  (for  the  moon),  k^  (for  the  sun),  or  (for 

both).  In  order  to  illustrate  the  error  due  to  fixing  e  in  the  model  in¬ 
volving  the  moon,  we  take  k  from  (5.2a)  and  consider  <j>  =  +%rr  resulting  in 
the  largest  magnitude  of  <N>  .  The  following  values  <N>  are  obtained  for 
chosen  e's: 

18.5°  ... -15.1  cm  ,  23.5°  ...  -13.6cm  ,  28.5°  ... -11.7 cm. 

Thus  the  worst  error  associated  with  e*23.5°  is  1.9cm;  for  <(>=0  this  er¬ 
ror  would  be  0.9cm.  The  acceptability  of  our  model  has  thus  been  establish¬ 
ed. 

The  relation  (5.9)  can  also  be  presented  as 

<N>  =  (3/8)k(cos  2<p- l/3)(cos  2e  +  1/3)  .  (5.10) 

This  formula  indicates  that  if  e  were  54.7°  (as  z'  in  equation  5.3)  the 
time  average  of  N  would  be  identically  zero  and  the  mean  sea  surface  wpald 
coincide  with  the  static  geoid.  If  the  x-axis  in  Figure  5.1  represented 
the  equator,  the  "ecliptic"  would  pass  through  the  point  P'.  As  a  matter 
of  interest  we  confirm  that  the  average  of  <N>  taken  over  the  sphere  is 
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zero,  due  to  cos  24>  =  1/3 ;  since  In  (5.5a)  we  had  N=0  at  any  given 
mement,  It  thus  follows  that 

<N>  =  <N>  =  0 

as  It  should  (the  order  of  Integration  does  not  matter). 

Upon  taking  e  =  23.5°  In  (5.10),  we  have 

<N>  =  k(cos  2<P-  1/3 )  ,  (5.11a) 

k  *  0.381k;  (5.11b) 

if  k  represents  kj  ,  k2  ,  or  kj  +  k2  ,  the  values  of  ic  are  0.102m,  0.047m, 
or  0.149m,  respectively.  In  terms  of  the  polar  distance  0( 9  =  -  <{>) , 

equation  (5.11a)  is  rewritten  as 

<N>  =  (-k)(cos  20  +  1/3)  . 

Thus  the  mean  sea  surface  in  our  model  is  a  geocentric  ellipsoid  of  revolu¬ 
tion  whose  minor  axis  passes  through  the  poles.  This  is  seen  immediately 
from  the  analogy  with  the  formula  (5.1)  and  Figure  5.1;  the  x-axis  and  the 
angle  z  are  now  replaced  by  the  polar  axis  and  the  angle  0,  and  the  negative 
value  of  (-k)  indicates  that  this  axis  is  the  minor  axis  of  an  ellipse. 

Upon  including  the  effect  of  both  the  moon  and  the  sun,  the  for¬ 
mulas  (5.11)  giving  the  separation  between  the  static  geoid  and  the  mean  sea 
surface  under  the  assumption  of  equilibrium  tide  is 

<N>  =  0.15x(cos  2<j>-l/3)  meter  .  (5.12) 
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From  (5.12)  the  largest  separation  is  found  for  the  poles  where  it  reaches 
-20cm.  For  the  equator  we  have  <N>  =  10cm.  This  or  a  similar  type  of  "cor¬ 
rection"  is  significant  if  a  consistent  relationship  between  the  two  sur¬ 
faces  is  needed  to  within  a  decimeter  precision. 
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5.3  Considerations  Related  to  Sea  Surface  Slopes 

It  has  been  mentioned  that  in  various  areas  the  sea  surface  may 
depart  from  an  equipotential  surface.  In  a  global  adjustment  of  satellite 
altimetry  the  latter  is  described  by  spherical  harmonic  potential  coeffici¬ 
ents  (C's  and  S's)  truncated  at  a  predetermined  degree  and  order;  for  ex¬ 
ample,  the  model  may  be  complete  through  the  degree  and  order  (14,14),  etc. 

If,  in  a  given  area,  the  sea  surface  is  believed  to  intersect  an  equipoten¬ 
tial  surface  along  a  line  of  known  direction,  a  parameter  describing  an 
average  slope  of  the  sea  surface  with  respect  to  the  equipotential  surface 
can  be  added  to  the  adjustable  C's  and  S's. 

There  may  be  as  many  new  parameters  of  this  type  as  there  are  areas 
suspect  to  exhibit  sea  surface  slopes.  One  such  area  can  be,  for  example, 
the  Equatorial  Pacific  Ocean  where  an  east-west  slope  may  be  indicated.  The 
line  of  intersection,  henceforth  called  "zero  line",  is  in  this  case  chosen 
to  run  north-south,  typically  through  the  middle  of  the  area  of  interest  (the 
area  can  be  delimited  by  the  geographic  coordinates  of  its  four  corners  and 
every  observation  point  is  checked  whether  it  falls  within  the  boundaries). 

However,  an  additional  parameter  may  be  included  in  an  adjustment, 
representing  a  vertical  shift  of  the  sea  surface.  This  has  the  effect  of  al¬ 
lowing  the  zero  line  to  translate,  parallel  to  itself,  to  the  right  or  left 
from  its  stipulated  position.  Accordingly,  two  sea  slope  parameters  (one  for 
the  slope  as  such,  one  for  the  shift)  per  region  of  interest  may  be  introduced 
in  a  global  adjustment. 
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The  slope  parameters  have  characteristics  of  an  average  with 
regard  to  both  the  surface  and  the  time.  The  first  property  follows  from 
the  fact  that  one  slope  parameter  applies  for  a  whole  area  and  indicates 
a  general  trend  of  the  sea  surface  from  one  end  of  the  area  to  the  other. 

And  the  second  property  can  be  ascribed  to  altimeter  data  being  usually 
collected  over  extended  periods  of  time  (e.g.,  several  months)  so  that  many 
periodic  (e.g.,  tides)  or  random  fluctuations  are  smoothed  out  by  the  ad¬ 
justment. 

If  a  time  variation  of  a  given  sea  slope  is  of  interest,  the  same 
kind  of  adjustment  may  be  repeated  with  data  sets  having  a  comparable spacial 
distribution  as  well  as  time  span,  but  collected  at  different  time  epochs. 

In  order  to  bring  the  sea  slope  parameters  into  focus  during  such  a  process, 
the  C's  and  S's  could  be  fixed  (eliminated  from  the  adjustment)  or  weighted 
at  the  values  which  would  be  the  same  during  all  such  adjustments. 

A  localized  adjustment  has  been  performed  in  the  past  with  point 
mass  magnitudes  as  parameters,  based  on  the  residuals  obtained  in  a  global 
adjustment  (via  C's  and  S's).  The  localized  adjustment  can  contain  the  sea 
slope  parameters  as  well,  but  there  is  likely  to  be  only  one  "area  of  in¬ 
terest",  i.e.,  the  total  area  of  the  point  mass  adjustment.  In  other  re¬ 
spects  the  sea  slope  parameters  can  be  treated  in  a  similar  way  as  described 
above.  Thus  the  localized  adjustment  may  contain  only  two  additional  para¬ 
meters,  provided  the  slope  and  the  (vertical)  shift  of  the  sea  surface  are 
indeed  sought  in  that  area. 


We  now  return  to  the  global  adjustment  and  illustrate  how  the  sea 
slope  parameters  could  be  useful  and  how  they  could  be  handled.  It  has  al¬ 
ready  been  suggested  that  they  may  contribute  to  a  "better"  adjustment.  In 
particular,  they  could  prevent  the  C's  and  S's  from  being  contaminated  by 
long-wavelength  distortions  represented  by  actual  sea  slopes  if  these  are 
significant.  Furthermore,  the  residuals  will  have  in  general  smaller  magni¬ 
tude  when  these  parameters  are  included  since  the  adjusted  surface  will  fit 
better  the  observed  surface  in  a  physically  meaningful  manner.  As  a  by¬ 
product,  know!  edge  about  sea  surface  slopes  can  serve,  in  its  own  right,  in 
specific  applications. 

An  example  of  two  areas  with  separate  sea  slope  parameters  is  il¬ 
lustrated  in  Figure  5.2.  An  added  complication  is  that  the  zero  line  is  not 
straight  which  accounts  for  an  auxiliary  area  wedged  between  the  other  two. 

In  this  particular  area  the  zero  line  reduces  to  a  point  (Q).  The  sea  slope 
parameters  in  this  area  could  be  taken  as  the  mean  of  the  corresponding  para¬ 
meters  in  the  neighboring  areas.  A  simpler  case  would  be  represented  by  a 
straight  zero  line  with  the  two  areas  separated  by  a  straight  line  (perpen¬ 
dicular  to  the  zero  line).  Eventually,  there  could  be  an  intermediate  area 
between  the  two  in  the  form  of  a  strip  (not  a  wedge)  where  the  values  of  the 
sea  slope  parameters  would  again  be  the  appropriate  means.  The  simplest  case 
consists  of  individual  disjoint  areas  or  of  only  one  area  attributed  sea 
slope  parameters. 


Illustration  of  two  basic  areas  (1  and  2)  and  a  bent  zero  line; 
the  area  1  is  delimited  by  the  points  1,  2,  3,  4,  Q 


The  adjusted  global  geoid  described  by  the  C's  and  S's  will  be 
smooth,  but  slight  discontinuities  in  the  residuals  may  be  experienced  at 
the  edges  of  a  sea  slope  area,  especially  far  from  the  zero  line,  except 
where  the  area  is  delimited  by  continents.  However,  this  represents  no 
cause  for  concern  since  a  subsequent  local  adjustment  usually  involves  a 
region  which  is  all  within  one  such  area,  and  even  if  there  should  be  an 
exception  to  this  statement  the  point  mass  parameters  would  provide  for  a 
smoothed  local  geoid. 
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If  the  sea  slope  areas  present  In  an  adjustment  cover  most  of 
the  world's  oceans  and  if  all  the  zero  lines  are  chosen  to  have  central 
locations,  one  would  expect  that 
all  areas 

I  (area),  x  (shift  parameter).  *  0  .  (5.13) 

i=l 

It  is  so  because  the  sea  surface  topography  is  not  likely  to  enclose  dif¬ 
ferent  volume  from  that  enclosed  by  an  equipotential  surface  (fluctuations 
of  the  sea  surface  around  an  equipotential  surface  are  believed  to  be  small 
and  to  have  random  characteristics).  In  fact,  equation  (5.13)  could  be 
adopted  as  a  constraint  or,  which  amounts  to  the  same  thing  in  practice,  as 
a  heavily  weighted  observation  equation. 

Before  forming  an  observation  equation  Involving  the  two  sea  slope 
parameters  ("t"  for  the  slope  Itself,  "q”  for  the  shift)  the  distance  between 
an  observation  point  and  the  zero  line  has  to  be  computed.  This  distance 
in  Figure  5.2)  will  have  the  sign  associated  with  It  in  such  a  way  that 
the  distance  is  positive  due  west  of  the  zero  line.  This  means  that  for  t>0, 
the  sea  surface  slopes  "up"  due  west  of  the  zero  line  and  "down"  due  east  of 
this  line;  for  t<0  the  situation  is  reversed.  Since,  In  the  observation 
equation,  the  constant  terms  are  expressed  in  meters  the  parameter  q  will 
have  units  m  and  the  parameter  t  will  have  units  imi/km  (i.e.,  meter/mega¬ 
meter),  provided  the  angular  distance  (in  radians)  is  multiplied  by  R,  the 
earth's  radius,  expressed  In  megameters  (R*  6.371  Mm).  The  value  of  l ^  is 
found  with  the  aid  of  Figure  5.3.  From  the  spherical  triangle  PAE  we  have 
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(5.14a) 


r 


sin  p  slni  =  sin  (Aq-A)  cose|)  , 
sin  p  cost  =  s1n<)>  ,  (5.14b) 

where  the  index  "o"  is  associated  with  the  zero  line  and  the  geographic 
coordinates  (♦♦A)  associated  with  the  observation  point  (A)  are  not  in¬ 
dexed.  From  AEA'  we  deduce 

sin  =  cosa  sin  p  cost  +  sinaQ  sin  p  sinT  .  (5.14c) 

When  (5.14a)  and  (5.14b)  are  used  in  (5.14c),  with  the  appropriate  sign 
is  given  as 

l A  =  arc  sin  [cosa0  sin<}>  +  sinaQ  sin  (Aq-A)  cos<j> ]  .  (5.15) 

The  values  of  cos  aQ  and  sinctQ  are  computed  once  and  for  all  in  the  region 
of  interest  (aQ  is  measured  as  indicated  in  Figure  5.3). 

In  case  of  the  situation  leading  to  £g  of  Figure  5.2,  a  simple 
application  of  the  law  of  cosine  yields 

cos£g  =  sind>  sin  <f>g  +  cos<b  cos  <j>q  cos  (Aq-a)  , 

sign  »  sign  (Ag  -  A)  , 

where  sin$p  and  cos  are  computed  once  and  for  all  for  the  whole  "wedge" 
region. 
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p 


Figure  5.3 

Spherical  triangles  depicting  the  zero  line  and  an  observation  point  (A) 


In  a  matrix  form,  the  "undulation"  in  one  area  due  to  the  sea 
slope  parameters  is 


AN  =  [1  IK]  pj  > 

where  l  can  be  of  the  type  l ^  or  l g  .  In  general,  this  can  be  written  as 
AN  =  [a  b]  j'qj  >  (5.16) 
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where 


in  the  area  of  interest 


(5.17a) 


a  =  0 
b  =  0 


outside  the  area  of  interest. 


(5.17b) 


If  there  are  more  than  one  area  of  interest,  they  are  numbered  1,2,...  The 
same  applies  for  the  coefficients  a,b  and  for  the  parameters  q,t  so  that  we 
may  have  aj  ,  b^  ;  a^  ,  *  •  •  •  Qj  *  tj  ;  ^  ,  tg  ;  •  • .  In  forming  the 

observation  equation  we  keep  in  mind  that  the  initial  values  of  the  para¬ 
meters  q,t  and  thus  of  AN  are  zero  so  that  the  constant  term  is  unchanged 
from  an  earlier  formulation  without  the  sea  surface  slopes.  The  part  cor¬ 
responding  to  (5.16)  simply  augments  the  "terrestrial  part"  in  the  original 
equation  of  satellite  altimetry  as  presented  in  Section  2.4  of  [Blaha,  1977]. 
With  the  other  symbols  described  —  and  expressed  in  detail  —  in  this  ref¬ 
erence,  we  present  schematically  the  augmented  observation  equation: 


(5.18) 


where 


al  ^1  *  a2  ^2  *  *** 


drc 

dC 

dC, 

dS 

<1 

*1 

*>2 


nO 


nm 


nm 


+  [(R°-r°  +  d)  -  Hl 


bj  are  determined  as  in  equations  (5.17). 
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APPENDIX 


COMPARISON  OF  THE  NODE-POINT  APPROACH 
AND  THE  POINT-MASS  APPROACH 


The  present  topic  deals  with  the  geoidal  adjustment  by  means  of 
node  point  parameters  and  point  mass  parameters.  The  node  point  approach 
has  been  treated  in  the  papers  [Hadgigeorge  and  Trotter,  1976],  abbreviated 
as  [HT],  and  [Eckhardt  and  Hadgigeorge,  1977],  abbreviated  as  [EH].  The 
notion  of  covariance  function  computed  via  the  spherical  harmonic  potential 
coefficients  is  used  in  both  papers;  however,  the  solution  bandwidth  in 
[HT]  corresponds  essentially  to  the  degrees  3  through  N  and  in  [EH] ,  to  the 
degrees  n  through  N,  where  n  was  chosen  as  7.  According  to  n,  the  reference 
surface  can  be  either  a  spheroid  (n=  3)  or  a  higher  order  surface  (n=7,etc.); 
in  general,  the  model  is  assumed  to  be  known  through  the  degree  n-  1.  The 
paper  [HT]  also  contains  the  spheroidal  multiquadric  surface  model  in  which 
the  covariance  function  is  replaced  by  the  "kernel  function". 

The  spacing  of  the  nodes  may  be  limited  by  the  data  distribution; 
in  any  event,  it  is  reasonable  to  require  that  at  least  two  observations  in 
each  direction  should  be  present  per  node  point.  The  node  points  are  assum¬ 
ed  to  be  fairly  uniformly  distributed  through  the  region  of  interest,  the 
nodal  spacing  being  s°  (degrees  of  arc).  If  the  solution  is  carried  out  in 
terms  of  spherical  harmonics,  the  harmonics  up  to  and  including  the  degree 
N=  180°/s°  could  be  adjusted  provided  the  global  distribution  of  data  exists 
and,  of  course,  a  sufficiently  large  computer  is  available.  The  covariance 
function  for  geoid  undulations  corresponding  to  the  bandwidth  n  through  N 


■  iii  »  i  'i-* 
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is  given  as 


N 

D<*  )  =  o2n  {N,N}  Pn(cos  f^)  ,  (A.  1) 

n=n 

where  D(^.j.)  represents  the  covariance  of  the  geoid  undulation  (denoted  N  — 
not  to  be  confused  with  "N"  of  the  truncation)  for  two  points,  i  and  j,  sep¬ 
arated  by  the  spherical  distance  ^ j ,  a *  {N  ,N>  is  the  degree  variance 
(degree  n)  for  geoid  undulations,  and  Pn(cos  ^j)  are  the  Legendre  functions 
in  argument  cos  The  cross-covariance  function  for  gravity  anomalies 

■  J 

(Ag)  and  geoid  undulations  (N)  can  similarly  be  written  as 

N 

H(^ij)  =  l  a*  (Ag.N)  Pn(cos  <1^.)  ,  (A. 2) 

n=n 

where  o 2  {Ag,N}  is  the  degree  variance  linking  gravity  anomalies  and  geoid 
undulations. 

In  the  above  papers  the  covariance  function  serves  to  express  the 

geoid  "undulation"  (N.)  at  a  geoidal  point  i  as  a  function  of  parameters  c. 

* 

at  selected  nodes  j,  namely 

Ni  =  l  D(,f,ij)cj  ’  (A'3) 

0 

where  the  summation  is  carried  out,  in  theory,  over  all  of  the  node  points; 

the  "undulation"  is  written  in  quotes  to  indicate  that  the  value  I'L  may  be 

referred  to  a  spheroid  (as  is  customary  for  geoid  undulations)  as  well  as 
to  a  higher  order  surface.  The  vector  of  parameters  c^  is  interpreted  as 

[cl  -  [MfNNTir'[Nl  ,  (A. 4) 
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where  [N]  is  the  vector  of  geoid  "undulations"  at  the  node  points  and  the 
matrix  [M{NNT }]  is  composed  of  variances  and  covariances  at  these  points 
computed  from  the  (single)  covariance  function,  such  as  (A.l)  in  the  present 
case.  This  follows  from  the  definition  of  covariance  function,  i.e.,  from 

D(*-j)  =  M{N.Nj}  , 

where  M  indicates  a  mean  value  over  a  given  domain  (in  the  present  context  the 
domain  for  the  covariance  functions  is  the  whole  globe  represented  by  a  unit 
sphere;  in  the  averaging  process  the  "unwanted"  harmonics  are  neglected). 

With  the  interpretation  (A. 4)  the  model  expressed  in  (A. 3)  is  the  classical 
prediction  model  where  the  values  N.  are  predicted  from  the  values  N.  the 
"best"  way,  in  the  sense  that  the  variances  of  the  prediction  computed  through 
the  same  operator  M  are  the  smallest  possible. 

In  the  following  step,  the  prediction  model  (A. 3)  is  adopted  as  the 
mathematical  model  to  be  used  in  a  least  squares  adjustment.  The  values  N.  , 
whose  number  is  much  higher  than  that  of  N-  ,  are  considered  as  observations 

J 

and  the  quantities  c.  are  adjusted  in  the  usual  sense  (VTPV =  minimum) .  In 

J 

this  process  D(ij>. .)  simply  represent  the  deterministic  partial  derivatives 

"  J 

and  c.  are  the  unknown  parameters.  These  parameters  can  then  be  used  to  esti- 

J 

mate  other  quantities,  not  only  the  geoid  "undulations"  at  arbitrary  points, 
but  also  gravity  "anomalies"  (the  same  meaning  is  attached  to  the  quotes  as 
previously),  etc.  The  formula  for  predicting  the  "anomalies"  is  again  based 
on  a  prediction  model  which,  in  this  case,  has  the  form 

Ag1  =  l  H(^j)cj  .  (A. 5) 
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When  considering  the  "kernel  function"  of  Section  2.2  of  [HT] , 
one  should  keep  in  mind  that  only  the  geoid  "undulations"  can  be  estimated 
from  the  adjustment  and  that  the  main  difference  between  that  approach  and 
the  one  with  the  covariance  function  consists  in  a  different  shape  of  the 
pertinent  functions  (kernel  versus  covariance  function).  Accordingly, 
only  the  covariance  function  approach  will  be  considered  in  the  upcoming 
comparison  with  the  point  mass  approach. 

The  mathematical  models  relating  the  geoid  "undulations"  and  gra¬ 
vity  "anomalies"  to  the  point  masses  can  be  obtained  from  (4.17)  and  (4.19) 
as 

Ni  =  ?  (c/^jHPMjj  ,  (A.6) 

J 

where  c  =  40.68x10®  m2 ,  is  the  separation  between  the  i-th  "undulation" 

*  J 

point  and  the  j-th  point  mass  and  (PM),  is  the  magnitude  of  the  j-th  point 

J 

mass  in  units  of  10"  the  earth's  mass;  and 

Ag.  =  l  (ab)ij(PM)j  ,  (A. 7) 

J 

where 

aij  =  C’/£ij  * 

c'  =  6.2565 xlO®  mgal  m  , 

bu  ■  Rbij«f  j  - 2  • 

bjj  =  R  -  R,  cos  ^  . 

R  =  earth's  mean  radius  (6.371x10®  m)  , 
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I<(  !<  «l  . 

d  =  depth  of  point  masses. 

If  properly  scaled,  the  model  (A. 6)  is  comparable  to  the  model 
(A. 3)  except  for  the  form  of  the  "covariance"  functions;  the  quotes  will 
be  used  to  include  the  function  c/t..  in  (A. 6)  which  has  not  been  derived 
as  a  covariance  function.  The  scaling  of  the  adjustment  models  (A. 3)  and 
(A. 6)  would  not  affect  the  respective  least  squares  solutions;  however, 
the  predictions  (A. 5)  and  (A. 7)  would  be  affected  if  the  scale  factors  did 
not  correspond  to  their  counterparts  in  the  original  models  (equations  A. 3 
and  A. 6).  In  the  sequel,  the  scaling  of  (A. 5)  and  (A. 7)  will  not  correspond 
to  that  of  (A. 3)  and  (A. 6)  for  practical  reasons;  however,  the  differences 
in  scaling  will  be  listed.  In  all  the  cases  considered,  the  scaling  will  be 
such  that  the  "covariance"  functions  in  (A. 3),  (A. 6)  as  well  as  in  (A. 5), 

(A. 7)  will  have  the  value  "100"  for  .  =  0.  Thus,  the  shape  of  the  functions 

*  J 

may  be  easily  compared  upon  first  inspection.  The  units  in  any  of  the  models 
present  no  cause  for  concern;  the  functions  may  be  imagined  scaled  also  in 
this  respect  so  that  the  units  are  comparable. 

The  cases  compared  correspond  to  the  spacing  of  node  points  in 
s°=6°  (equilateral)  intervals  which  are  considered  to  entail  the  degree  and 
order  truncation  N=30.  The  covariance  functions  for  this  N  are  then  com¬ 
puted  according  to  (A.l)  and  (A. 2),  for  the  cases  n=3  and  n=7.  The  hori¬ 
zontal  separation  of  point  masses  is  also  6°  or  .667x10s  m.  For  illustra¬ 
tion  purposes  three  different  depths  of  the  point  masses  have  been  selected: 
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d  =  1.6  s  =  1.07X10*  m  , 

(A. 8a) 

d  =  .52s  =  0.35x10*  m  , 

(A. 8b) 

d  =  .80s  =  0.53x10*  m  . 

(A. 8c) 

Originally,  the  ratio  d/s  =  0.8  as  in  (A. 8c)  has  been  suggested  in  [Blaha, 
1977]  under  the  assumption  that  the  spacing  s°  of  point  masses  should  pro¬ 
vide  an  adequate  description  of  the  geoid  to  within  the  degree  and  order 
N=180°/s°.  Later,  in  [Blaha,  1979],  a  one-decimeter  precision  of  the  ad¬ 
justed  geoid  was  contemplated;  a  recommendation  was  reached  to  the  effect 
that  the  number  of  point  masses  should  be  doubled  in  each  direction  while 
their  depth  should  remain  unchanged.  This  would  result  in  the  angular  sep¬ 
aration  of  point  masses  s°  =  %180°/N  and  in  the  ratio  d/s=  1.6.  However,  to 
make  a  meaningful  comparison  with  the  present  node  point  approach,  the  orig¬ 
inal  arrangement  has  been  reverted  to  in  (A. 8c);  in  particular,  s°=180°/30 
*  6°  has  been  adhered  to  as  well  as  d/s *0.8.  If  a  more  faithful  representa¬ 
tion  of  the  geoid  were  required,  the  number  of  point  masses  as  well  as  the 
number  of  node  points  would  be  doubled  along  each  dimension.  But  in  that 
case  the  comparisons  would  be  completely  analogous  to  those  carried  out  pre¬ 
sently. 

The  "covariance"  functions  corresponding  to  the  adjustment  models 
(A. 3)  and  (A. 6)  are  presented  in  Figure  A.l.  The  cases  denoted  1)  and  2) 
represent  the  (scaled)  covariance  function  of  the  node  point  approach;  the 
curve  1)  is  the  result  of  the  choice  n  =  3  while  the  curve  2)  is  the  result 
of  the  choice  n=  7.  For  higher  and  higher  order  reference  surface  (i.e., 
for  greater  values  of  n)  the  covariance  function  becomes  markedly  steeper. 


The  cases  3),  4),  5)  in  the  figure  represent  the  "covariance"  function  of 
the  point  mass  approach  and  they  correspond  to  the  depths  of  point  masses 
selected  in  (A. 8a),  (A. 8b)  and  (A. 8c),  respectively.  The  shallower  the 
point  masses,  the  steeper  the  "covariance"  function.  In  either  approach, 
when  this  function  is  steeper  the  local  detail  can  be  better  represented 
provided  the  data  are  sufficient. 

The  "cross-covariance"  functions  corresponding  to  the  prediction 
models  (A. 5)  and  (A. 7)  appear  in  Figure  A. 2  where  the  same  five  cases  are 
presented.  All  the  curves  are  now  steeper  than  their  counterparts  in  Fig¬ 
ure  A.l;  this  is  natural  since  gravity  anomalies  represent  a  less  smooth 
function  than  geoid  undulations  (the  same  holds  true  for  "anomalies"  versus 
"undulations").  As  in  Figure  A.l,  greater  values  of  n  and  smaller  values 
of  d  give  rise  to  steeper  functions.  In  both  figures,  < p°  is  the  angle 
^ j  in  degrees.  Below  are  listed  the  numerical  scale  factors  used  for  the 
geoid  "undulation"  models  (Figure  A.l)  and  the  gravity  "anomaly"  models 
(Figure  A. 2): 

Case  1)  Scale  factor  6,042  (forN)  and  3.397  (forAg),  ratio  0.5622, 

Case  2)  Scale  factor  0.837  (forN)  and  1.336  (forAg),  ratio  1.596; 

Case  3)  Scale  factor  0.3811  (forN)  and  0.2325  (for  Ag) ,  ratio  0.6101, 

Case  5)  Scale  factor  0.7621  (forN)  and  1.164  (forAg),  ratio  1.528, 

Case  4)  Scale  factor  1.162  (forN)  and  2.896  (forAg),  ratio  2.492. 

The  above  numbering  of  cases  3,  5  and  4  corresponds  to  their  degree  of  steep¬ 
ness.  When  the  above  scale  factors  are  applied  to  the  functions  of  the  fig¬ 
ures,  original  (unsealed)  numerical  values  of  the  functions  are  recovered. 
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The  term  "ratio"  indicates  the  ratio  of  the  two  scale  factors  in  each  case 
(scale  factor  for  Ag  divided  by  the  scale  factor  for  N).  In  both  approaches 
a  steeper  curve  exhibits  a  higher  "ratio". 

Initially,  the  case  5)  of  the  point  mass  approach  was  to  be  com¬ 
pared  with  the  cases  1)  and  2)  of  the  node  point  approach  corresponding 
respectively  to  the  adjustments  carried  out  in  [HT]  and  [EH].  From  the  fig¬ 
ures  it  appears  that  the  case  5)  is  quite  comparable  with  the  case  2)  where 
n=  7.  For  the  sake  of  a  better  understanding  of  the  point  mass  approach 
two  additional  cases  have  been  included,  one  corresponding  to  deeper  masses 
(case  3)  and  the  other  corresponding  to  shallower  point  masses  (case  4). 

These  comparisons  have  indicated  that  the  main  difference  between 
the  node-point  approach  and  the  point-mass  approach  consists  in  somewhat  dif¬ 
ferent  shapes  of  the  "covariance"  functions.  In  certain  cases  the  disparity 
of  the  two  approaches  is  not  great;  this  is  apparent  for  the  cases  2)  ver¬ 
sus  5)  where  a  fairly  good  agreement  may  be  observed  in  the  two  figures, 
especially  for  smaller  values  of  t|>°.  Here  the  covariance  function  construct¬ 
ed  for  the  degrees  n  =  7  and  N*  30  corresponds  to  the  depth  0.53X106  m  of  the 
points  masses  distributed  exactly  as  the  node  points.  A  remarkable  agree¬ 
ment  is  noticed  with  respect  to  the  values  "ratio"  for  these  two  cases.  This 
seems  to  indicate  that  not  only  the  adjusted  values  (geoid  "undulations") 
but  also  the  predicted  values  (gravity  "anomalies"  in  addition  to  predicted 
"undulations"  at  arbitrary  locations)  could  agree  fairly  closely  between 
the  two  adjustments.  In  other  instances  the  difference  between  the  two  ap¬ 
proaches  may  become  significant,  for  example  if  shallower  point  masses  are 
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not  matched  by  a  sufficiently  high-order  reference  surface.  Comparisons 
with  other  localized  adjustment  models  would  proceed  along  similar  lines 
with  emphasis  on  the  corresponding  "covariance"  functions,  the  earlier- 
mentioned  multiquadric  model  being  a  good  example  in  this  respect. 
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